—
—

              
package ICC::Support::spline;
use strict;
use Carp;
our $VERSION = 0.12;
# revised 2018-09-04
#
# Copyright © 2004-2020 by William B. Birkett
# 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', 'output', 'type', 'damp', 'fix_hl', 'fix_sh'
# default values for 'input' and 'output' are [0, 1]
# parameters: ([ref_to_attribute_hash])
# returns: (ref_to_object)
sub new {
        # get object class
        my $class = shift();
        # create empty spline object
        my $self = [
                {}, # object header
                [], # input values
                [], # 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));
}
# create inverse object
# swaps the input and output arrays
# returns: (ref_to_object)
sub inv {
        # clone new object
        my $inv = Storable::dclone(shift());
        # get min/max values
        my @s = _minmax($inv);
        # verify object is monotonic
        (! @s) or croak("cannot invert a non-monotonic 'spline' object");
        # if input is a range
        if (@{$inv->[1]} == 2 && @{$inv->[2]} > 2) {
                 
                # interpolate input array
                $inv->[1] = [map {(1 - $_/$#{$inv->[2]}) * $inv->[1][0] + $_/$#{$inv->[2]} * $inv->[1][-1]} (0 .. $#{$inv->[2]})];
                 
        }
        # swap input and output arrays
        ($inv->[1], $inv->[2]) = ($inv->[2], $inv->[1]);
        # clear parameteric derivative matrix
        $inv->[5] = [];
        # sort input/output values
        _sort_values($inv);
        # compute knot derivatives
        _objderv($inv);
        # add segment min/max y-values
        _minmax($inv);
        # return inverted object
        return($inv);
}
# 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 {
        # verify 4 parameters
        (@_ == 4) or croak('wrong number of parameters');
        # write spline data to profile
        goto &_writeICCspline;
}
# 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) = @_;
        # local variable
        my ($low);
        # if an extrapolated solution (s < 0)
        if ($in < $self->[1][0]) {
                 
                # return extrapolated value
                return($self->[2][0] + $self->[3][0] * ($in - $self->[1][0]));
                 
        # if an extrapolated solution (s >= 1)
        } elsif ($in >= $self->[1][-1]) {
                 
                # return extrapolated value
                return($self->[2][-1] + $self->[3][-1] * ($in - $self->[1][-1]));
                 
        # if input array contains x-values
        } elsif ($#{$self->[1]} > 1) {
                 
                # initilize segment
                $low = 0;
                 
                # while input value > upper knot value
                while ($in > $self->[1][$low + 1]) {
                         
                        # increment segment
                        $low++;
                         
                }
                 
                # return interpolated value
                return(_fwd($self, ($self->[1][$low] - $in)/($self->[1][$low] - $self->[1][$low + 1]), $low));
                 
        # input array contains range
        } 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 ($ix, @sol, @t);
        # get output array upper index
        $ix = $#{$self->[2]};
        # 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] + ($in - $self->[2][0])/$self->[3][0]);
                 
        }
        # for each segment (0 <= s < 1)
        for my $i (0 .. $ix - 1) {
                 
                # if input is lower knot
                if ($in == $self->[2][$i]) {
                         
                        # add solution (x-values : range)
                        push(@sol, ($#{$self->[1]} > 1) ? $self->[1][$i] : ($self->[1][0] * ($ix - $i) + $self->[1][1] * $i)/$ix);
                         
                }
                 
                # 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);
                         
                        # if input array contains x-values
                        if ($#{$self->[1]} > 1) {
                                 
                                # add solution(s)
                                push(@sol, map {(1 - $_) * $self->[1][$i] + $_ * $self->[1][$i + 1]} @t);
                                 
                        # input array contains range
                        } else {
                                 
                                # add solution(s)
                                push(@sol, map {($self->[1][0] * ($ix - $i - $_) + $self->[1][1] * ($i + $_))/$ix} @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] + ($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) = @_;
        # local variable
        my ($low);
        # if an extrapolated solution (s < 0)
        if ($in < $self->[1][0]) {
                 
                # return endpoint derivative
                return($self->[3][0]);
                 
        # if an extrapolated solution (s >= 1)
        } elsif ($in >= $self->[1][-1]) {
                 
                # return endpoint derivative
                return($self->[3][-1]);
                 
        # if input array contains x-values
        } elsif ($#{$self->[1]} > 1) {
                 
                # initilize segment
                $low = 0;
                 
                # while input value > upper knot value
                while ($in > $self->[1][$low + 1]) {
                         
                        # increment segment
                        $low++;
                         
                }
                 
                # return interpolated derivative value
                return(_derv($self, ($self->[1][$low] - $in)/($self->[1][$low] - $self->[1][$low + 1]), $low));
                 
        # input array contains range
        } else {
                 
                # return interpolated derivative value
                return(_derv($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($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 ($dx, $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 ($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] * ($in - $self->[1][0]);
                         
                }
                 
        # if an extrapolated solution (s >= 1)
        } elsif ($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] * ($in - $self->[1][-1]);
                         
                }
                 
        } else {
                 
                # if input array contains x-values
                if ($#{$self->[1]} > 1) {
                         
                        # initialize segment
                        $low = 0;
                         
                        # while input value > upper knot value
                        while ($in > $self->[1][$low + 1]) {
                                 
                                # increment segment
                                $low++;
                                 
                        }
                         
                        # compute delta x-value
                        $dx = $self->[1][$low + 1] - $self->[1][$low];
                         
                        # compute t-value
                        $t = ($in - $self->[1][$low])/$dx;
                         
                # input array contains range
                } else {
                         
                        # compute delta x-value
                        $dx = ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
                         
                        # compute t-value, index
                        ($t, $low) = POSIX::modf(($in - $self->[1][0])/$dx);
                         
                }
                 
                # 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 * $dx * $self->[5][$low][$i] + $h11 * $dx * $self->[5][$low + 1][$i];
                         
                }
                 
        }
        # return
        return(@pj);
}
# get/set header hash reference
# 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 array reference
# an array of two values is a range
# otherwise, it contains input values
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub input {
        # get object reference
        my $self = shift();
        # if one parameter supplied
        if (@_ == 1) {
                 
                # verify 'input' vector (two or more numeric values)
                (ICC::Shared::is_num_vector($_[0]) && 2 <= @{$_[0]}) or croak('invalid input values');
                 
                # verify 'input' vector size
                (@{$_[0]} == 2 || @{$_[0]} == @{$self->[2]}) or carp('number of input values differs from object output values');
                 
                # store output values
                $self->[1] = [@{$_[0]}];
                 
                # update object internal elements
                update($self, 1);
                 
        } elsif (@_) {
                 
                # error
                carp('too many parameters');
                 
        }
        # return array reference
        return($self->[1]);
}
# get/set output array reference
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub output {
        # get object reference
        my $self = shift();
        # if one parameter supplied
        if (@_ == 1) {
                 
                # verify 'output' vector (two or more numeric values)
                (ICC::Shared::is_num_vector($_[0]) && 2 <= @{$_[0]}) or croak('invalid output values');
                 
                # verify 'output' vector size
                (@{$self->[1]} == 2 || @{$_[0]} == @{$self->[1]}) or carp('number of output values differs from object input values');
                 
                # store output values
                $self->[2] = [@{$_[0]}];
                 
                # update object internal elements
                update($self);
                 
        } elsif (@_) {
                 
                # error
                carp('too many parameters');
                 
        }
        # return array reference
        return($self->[2]);
}
# normalize transform
# adjusts Bernstein coefficients linearly
# adjusts 'input' and 'output' by default
# hash keys: 'input', 'output'
# parameter: ([hash_ref])
# returns: (ref_to_object)
sub normalize {
        # get parameters
        my ($self, $hash) = @_;
        # local variables
        my ($m, $b, $val, $src, $x0, $x1, $y0, $y1, $f0, $f1, @minmax, $p);
        # set hash key array
        my @key = qw(input output);
        # for ('input', 'output')
        for my $i (0, 1) {
                 
                # undefine slope
                $m = undef;
                 
                # if no parameter (default)
                if (! defined($hash)) {
                         
                        # if array has 2 or more coefficients
                        if (@{$self->[$i + 1]} > 1 && ($x0 = $self->[$i + 1][0]) != ($x1 = $self->[$i + 1][-1])) {
                                 
                                # special case, if 'output' and 'fix_hl' was called
                                $x0 = 0 if ($i && defined($self->[0]{'hl_retention'}));
                                 
                                # compute adjustment values
                                $m = abs(1/($x0 - $x1));
                                $b = -$m * ($x0 < $x1 ? $x0 : $x1);
                                 
                                # set endpoints
                                ($y0, $y1) = $x0 < $x1 ? (0, 1) : (1, 0);
                                 
                        }
                         
                } else {
                         
                        # if hash contains key
                        if (defined($val = $hash->{$key[$i]})) {
                                 
                                # if array has 2 or more coefficients
                                if (@{$self->[$i + 1]} > 1 && $self->[$i + 1][0] != $self->[$i + 1][-1]) {
                                         
                                        # if hash value is an array reference
                                        if (ref($val) eq 'ARRAY') {
                                                 
                                                # if two parameters
                                                if (@{$val} == 2) {
                                                         
                                                        # get parameters
                                                        ($y0, $y1) = @{$val};
                                                         
                                                        # get endpoints
                                                        $x0 = $self->[$i + 1][0];
                                                        $x1 = $self->[$i + 1][-1];
                                                         
                                                # if three parameters
                                                } elsif (@{$val} == 3) {
                                                         
                                                        # get parameters
                                                        ($src, $y0, $y1) = @{$val};
                                                         
                                                        # if source is 'endpoints'
                                                        if (! ref($src) && $src eq 'endpoints') {
                                                                 
                                                                # get endpoints
                                                                $x0 = $self->[$i + 1][0];
                                                                $x1 = $self->[$i + 1][-1];
                                                                 
                                                        # if source is 'minmax'
                                                        } elsif (! ref($src) && $src eq 'minmax') {
                                                                 
                                                                # if 'output' curve
                                                                if ($i) {
                                                                         
                                                                        # get all min/max values
                                                                        @minmax = (map {@{$_}} @{$self->[4]});
                                                                         
                                                                        # get min/max values
                                                                        $x0 = List::Util::min(@minmax);
                                                                        $x1 = List::Util::max(@minmax);
                                                                         
                                                                # if 'input' curve
                                                                } else {
                                                                         
                                                                        # get min/max values (endpoints)
                                                                        $x0 = $self->[1][0];
                                                                        $x1 = $self->[1][-1];
                                                                         
                                                                }
                                                                 
                                                        }
                                         
                                                # if four parameters
                                                } elsif (@{$val} == 4) {
                                                         
                                                        # get parameters
                                                        ($x0, $x1, $y0, $y1) = @{$val};
                                                         
                                                }
                                                 
                                                # if x/y values are numeric
                                                if ((4 == grep {Scalar::Util::looks_like_number($_)} ($x0, $x1, $y0, $y1)) && ($x0 != $x1)) {
                                                         
                                                        # compute slope
                                                        $m = ($y0 - $y1)/($x0 - $x1);
                                                         
                                                        # compute offset
                                                        $b = $y0 - $m * $x0;
                                                         
                                                } else {
                                                         
                                                        # warning
                                                        carp("invalid '$key[$i]' parameter(s)");
                                                         
                                                }
                                         
                                        }
                                 
                                } else {
                                         
                                        # warning
                                        carp("can't normalize '$key[$i]' coefficients");
                                         
                                }
                                 
                        }
                 
                }
                 
                # if slope is defined
                if (defined($m)) {
                         
                        # set endpoint flags
                        $f0 = $self->[$i + 1][0] == $x0;
                        $f1 = $self->[$i + 1][-1] == $x1;
                         
                        # adjust Bernstein coefficients (y = mx + b)
                        @{$self->[$i + 1]} = map {$m * $_ + $b} @{$self->[$i + 1]};
                         
                        # set endpoints (to ensure exact values)
                        $self->[$i + 1][0] = $y0 if ($f0);
                        $self->[$i + 1][-1] = $y1 if ($f1);
                         
                        # save input polarity
                        $p = $m < 0 if ($i == 0);
                         
                        # if 'input' and 'hl_retention' is defined
                        if (! $i && defined($val = $self->[0]{'hl_retention'})) {
                                 
                                # adjust value
                                $self->[0]{'hl_retention'} = $m * $val + $b;
                                 
                        }
                         
                        # if 'output' and 'hl_original' is defined
                        if ($i && defined($val = $self->[0]{'hl_original'})) {
                                 
                                # adjust values
                                $self->[0]{'hl_original'} = [map {$m * $_ + $b} @{$val}];
                                 
                        }
                         
                        # if 'output' and 'sh_original' is defined
                        if ($i && defined($val = $self->[0]{'sh_original'})) {
                                 
                                # adjust values
                                $self->[0]{'sh_original'} = [map {$m * $_ + $b} @{$val}];
                                 
                        }
                         
                }
                 
        }
        # update object internals, sorting values if 'input' polarity changed
        update($self, $p);
        # return
        return($self);
}
# check if monotonic
# returns min/max values
# returns s-values by default
# returns input-values if format is 'input'
# returns output-values if format is 'output'
# curve is monotonic if no min/max values
# parameters: ([format])
# returns: (array_of_values)
sub monotonic {
        # get object reference
        my ($self, $fmt) = @_;
         
        # local variables
        my ($t, $low);
        # if format undefined
        if (! defined($fmt)) {
                 
                # return s-values
                return(_minmax($self));
                 
        # if format 'input'
        } elsif ($fmt eq 'input') {
                 
                # if input array contains x-values
                if ($#{$self->[1]} > 1) {
                         
                        # return input values
                        return(map {($t, $low) = POSIX::modf($_); (1 - $t) * $self->[1][$low] + $t * $self->[1][$low + 1]} _minmax($self));
                         
                # input array contains range
                } else {
                         
                        # return input values
                        return(map {$t = $_/$#{$self->[2]}; (1 - $t) * $self->[1][0] + $t * $self->[1][1]} _minmax($self));
                         
                }
                 
        # if format 'output'
        } elsif ($fmt eq 'output') {
                 
                # return output values
                return(map {_fwd($self, POSIX::modf($_))} _minmax($self));
                 
        } else {
                 
                # error
                croak('unsupported format for min/max values');
                 
        }
         
}
# update internal object elements
# this method used primarily when optimizing
# call after changing 'input' or 'output' values
# if the 'input' values were changed, set the flag
# parameter: ([flag])
sub update {
        # get parameters
        my ($self, $flag) = @_;
        # sort values if flag set
        _sort_values($self) if ($flag);
        # recompute knot derivatives
        _objderv($self);
        # recompute segment min/max y-values
        _minmax($self);
        # recompute parametric derivative matrix if flag set, or 'akima' type spline
        _objpara($self) if ($flag || (defined($self->[0]{'type'} && $self->[0]{'type'} eq 'akima')));
}
# 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) or 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 equivalent '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(@_)));
}
# 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, $fmt2);
        # 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);
        # if input_parameter_array defined
        if (defined($self->[1])) {
                 
                # set format
                $fmt2 = join(', ', ('%.3f') x @{$self->[1]});
                 
                # append parameter string
                $s .= sprintf("  input parameters  [$fmt2]\n", @{$self->[1]});
                 
        }
        # if output_parameter_array defined
        if (defined($self->[2])) {
                 
                # set format
                $fmt2 = join(', ', ('%.3f') x @{$self->[2]});
                 
                # append parameter string
                $s .= sprintf("  output parameters [$fmt2]\n", @{$self->[2]});
                 
        }
        # return
        return($s);
}
# combined transform
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, 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));
                 
        }
         
}
# combined derivative
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, 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));
                 
        }
         
}
# compute knot derivatives
# parameter: (ref_to_object)
sub _objderv {
        # get parameter
        my $self = shift();
        # verify input values are unique
        _unique($self->[1]) or croak('input values not unique');
         
        # if spline type undefined or 'natural'
        if (! defined($self->[0]{'type'}) || $self->[0]{'type'} eq 'natural') {
                 
                # compute knot derivatives for natural spline
                _natural_derv($self)
                 
        # if spline type is 'akima'
        } elsif ($self->[0]{'type'} eq 'akima') {
                 
                # compute knot derivatives for Akima spline
                _akima_derv($self);
                 
        } else {
                 
                # error
                croak('undefined spline type');
                 
        }
         
}
# compute knot derivatives for natural spline
# parameter: (ref_to_object)
sub _natural_derv {
        # get parameter
        my $self = shift();
        # local variables
        my ($ix, $rhs, $dx, $info, $derv);
        # verify object has two or more knots
        (($ix = $#{$self->[2]}) > 0) or 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');
        # if input array contains input values
        if ($#{$self->[1]} > 1) {
                 
                # for each inner element
                for my $i (1 .. $ix - 1) {
                         
                        # compute rhs (6 * (y[i + 1] - y[i - 1])/(x[i + 1] - x[i - 1]))
                        $rhs->[$i][0] = 6 * ($self->[2][$i + 1] - $self->[2][$i - 1])/($self->[1][$i + 1] - $self->[1][$i - 1]);
                         
                }
                 
                # set rhs endpoint values
                $rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0])/($self->[1][1] - $self->[1][0]);
                $rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1])/($self->[1][$ix] - $self->[1][$ix - 1]);
                 
        # input array contains range
        } else {
                 
                # compute x-segment length
                $dx = ($self->[1][1] - $self->[1][0])/$ix;
                 
                # for each inner element
                for my $i (1 .. $ix - 1) {
                         
                        # compute rhs (6 * (y[i + 1] - y[i - 1]))/(x[i + 1] - x[i - 1])
                        $rhs->[$i][0] = 3 * ($self->[2][$i + 1] - $self->[2][$i - 1])/$dx;
                         
                }
                 
                # set rhs endpoint values
                $rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0])/$dx;
                $rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1])/$dx;
                 
        }
        # 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 knot derivatives for Akima spline
# parameter: (ref_to_object)
sub _akima_derv {
        # get parameter
        my $self = shift();
        # local variables
        my ($ix, @m, $s, $b, $d1, $d2);
        # verify object has two or more knots
        (($ix = $#{$self->[2]}) > 0) or croak('object must have two or more knots');
        # compute uniform segment length (for 'input' range)
        $s = ($self->[1][1] - $self->[1][0])/$ix if ($#{$self->[1]} == 1);
        # for each segment
        for my $i (0 .. ($ix - 1)) {
                 
                # compute segment slope (range // x-values)
                $m[$i + 2] = ($self->[2][$i + 1] - $self->[2][$i])/($s // ($self->[1][$i + 1] - $self->[1][$i]));
                 
        }
        # if 2 points
        if ($ix == 1) {
                 
                # linear slopes
                $m[0] = $m[1] = $m[3] = $m[4] = $m[2];
                 
        # if 3 or more points
        } else {
                 
                # quadratic slopes
                $m[1] = 2 * $m[2] - $m[3];
                $m[0] = 2 * $m[1] - $m[2];
                $m[$ix + 2] = 2 * $m[$ix + 1] - $m[$ix];
                $m[$ix + 3] = 2 * $m[$ix + 2] - $m[$ix + 1];
                 
        }
        # if damping value is undefined
        if (! defined($self->[0]{'damp'})) {
                 
                # set default damping value
                $self->[0]{'damp'} = List::Util::sum(map {$_ * $_} @m)/(@m * 50);
                 
        }
        # Akima used abs() functions to calulate the slopes
        # this can lead to abrupt changes in curve shape and parametric derivatives
        # we replace abs(x) with sqrt(x^2 + b) to add stability
        # get damping value and verify >= 0
        (($b = $self->[0]{'damp'}) >= 0) or carp('akima damping value is negative');
        # compute Akima derivative
        for my $i (0 .. $ix) {
                 
                # compute slope differences (with damping)
                $d1 = sqrt(($m[$i + 1] - $m[$i])**2 + $b);
                $d2 = sqrt(($m[$i + 3] - $m[$i + 2])**2 + $b);
                 
                # if denominator not 0
                if ($d1 + $d2) {
                         
                        # use Akima slope (weighted average with damping)
                        $self->[3][$i] = ($d2 * $m[$i + 1] + $d1 * $m[$i + 2])/($d1 + $d2);
                         
                } else {
                         
                        # otherwise, use average slope of adjoining segments
                        $self->[3][$i] = ($m[$i + 1] + $m[$i + 2])/2;
                         
                }
                 
        }
         
}
# compute partial derivative matrix
# parameter: (ref_to_object)
sub _objpara {
        # get parameter
        my $self = shift();
        # verify input values are unique
        _unique($self->[1]) or croak('input values not unique');
         
        # if spline type undefined or 'natural'
        if (! defined($self->[0]{'type'}) || $self->[0]{'type'} eq 'natural') {
                 
                # compute knot derivatives for natural spline
                _natural_para($self)
                 
        # if spline type is 'akima'
        } elsif ($self->[0]{'type'} eq 'akima') {
                 
                # compute knot derivatives for Akima spline
                _akima_para($self);
                 
        } else {
                 
                # error
                croak('undefined spline type');
                 
        }
         
}
# compute partial derivative matrix for natural spline
# parameter: (ref_to_object)
sub _natural_para {
        # get parameter
        my $self = shift();
        # local variables
        my ($ix, $n, $rhs, $x, $info, $derv);
        # verify object has two or more knots
        (($ix = $#{$self->[2]}) > 0) or 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');
        # if input array contains input values
        if ($#{$self->[1]} > 1) {
                 
                # for each inner element
                for my $i (1 .. $ix - 1) {
                         
                        # set rhs diagonal values
                        $rhs->[$i][$i - 1] = $x = 6/($self->[1][$i - 1] - $self->[1][$i + 1]);
                        $rhs->[$i][$i + 1] = -$x;
                         
                }
                 
                # set rhs endpoint values
                $rhs->[0][0] = $x = 3/($self->[1][0] - $self->[1][1]);
                $rhs->[0][1] = -$x;
                $rhs->[$ix][$ix] = $x = 3/($self->[1][$ix] - $self->[1][$ix - 1]);
                $rhs->[$ix][$ix -1] = -$x;
                 
        # input array contains range
        } else {
                 
                # compute rhs value
                $x = 3 * $ix/($self->[1][1] - $self->[1][0]);
                 
                # for each inner element
                for my $i (1 .. $ix - 1) {
                         
                        # set rhs diagonal values
                        $rhs->[$i - 1][$i] = $x;
                        $rhs->[$i + 1][$i] = -$x;
                         
                }
                 
                # set rhs endpoint values
                $rhs->[0][0] = -$x;
                $rhs->[1][0] = -$x;
                $rhs->[$ix][$ix] = $x;
                $rhs->[$ix - 1][$ix] = $x;
                 
        }
        # 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();
                 
        }
         
}
# compute partial derivative matrix for Akima spline
# parameter: (ref_to_object)
sub _akima_para {
        # get object
        my $self = shift();
        # local variables
        my ($ix, $s, @m, @x, $b);
        my ($d1, $d2, $dif, @dm, $dd1, $dd2, $derv);
         
        # see 'akima_parametric.plx' for details on this function
        # verify object has two or more knots
        (($ix = $#{$self->[2]}) > 0) or croak('object must have two or more knots');
        # compute uniform segment length (for 'input' range)
        $s = ($self->[1][1] - $self->[1][0])/$ix if ($#{$self->[1]} == 1);
        # for each segment
        for my $i (0 .. ($ix - 1)) {
                 
                # compute segment x-length
                $x[$i + 2] = $s // ($self->[1][$i + 1] - $self->[1][$i]);
                 
                # compute segment slope (range // x-values)
                $m[$i + 2] = ($self->[2][$i + 1] - $self->[2][$i])/$x[$i + 2];
                 
        }
        # if 2 points
        if ($ix == 1) {
                 
                # linear slopes
                $m[0] = $m[1] = $m[3] = $m[4] = $m[2];
                 
        # if 3 or more points
        } else {
                 
                # quadratic slopes
                $m[1] = 2 * $m[2] - $m[3];
                $m[0] = 2 * $m[1] - $m[2];
                $m[$ix + 2] = 2 * $m[$ix + 1] - $m[$ix];
                $m[$ix + 3] = 2 * $m[$ix + 2] - $m[$ix + 1];
                 
        }
        # get the damping value
        $b = $self->[0]{'damp'};
        # for each node
        for my $i (0 .. $ix) {
                 
                # add row of zeros
                $derv->[$i] = [(0) x ($ix + 1)];
                 
        }
        # if 2 knots
        if ($ix == 1) {
                 
                # set partial derivatives
                $derv = [[-1/$x[2], 1/$x[2]], [-1/$x[2], 1/$x[2]]];
                 
        } else {
                 
                # for each row
                for my $i (0 .. $ix) {
                         
                        # compute slope differences (with damping)
                        $d1 = sqrt(($m[$i + 1] - $m[$i])**2 + $b);
                        $d2 = sqrt(($m[$i + 3] - $m[$i + 2])**2 + $b);
                         
                        # for each column
                        for my $j (0 .. $ix) {
                                 
                                # compute difference
                                $dif = $j - $i;
                                 
                                # skip outer corner regions
                                next if (abs($dif) > 2);
                                 
                                # the @dm array contains (∂m0/∂n, ∂m1/∂n, ∂m2/∂n, ∂m3/∂n)
                                # where m0 - m3 are the chord slopes, and n is the node value
                                 
                                # if difference -2
                                if ($dif == -2) {
                                         
                                        # if row 1
                                        if ($i == $ix) {
                                                 
                                                # set partial derivatives
                                                @dm = (-1/$x[$i], 0, 1/$x[$i], 2/$x[$i]);
                                                 
                                        } else {
                                                 
                                                # set partial derivatives
                                                @dm = (-1/$x[$i], 0, 0, 0);
                                                 
                                        }
                                         
                                # if difference -1
                                } elsif ($dif == -1) {
                                         
                                        # if row 1
                                        if ($i == 1) {
                                                 
                                                # if more than 3 knots
                                                if ($ix > 2) {
                                                         
                                                        # set partial derivatives
                                                        @dm = (-2/$x[$i + 1], -1/$x[$i + 1], 0, 0);
                                                         
                                                } else {
                                                         
                                                        # set partial derivatives
                                                        @dm = (-2/$x[$i + 1], -1/$x[$i + 1], 0, 1/$x[$i + 1]);
                                                         
                                                }
                                                 
                                        # if next to last row
                                        } elsif ($i == $ix - 1) {
                                                 
                                                # set partial derivatives
                                                @dm = (1/$x[$i], -1/$x[$i + 1], 0, 1/$x[$i + 1]);
                                                 
                                        # if last row
                                        } elsif ($i == $ix) {
                                                 
                                                # set partial derivatives
                                                @dm = (1/$x[$i], -1/$x[$i + 1], -2/$x[$i + 1] - 1/$x[$i], -3/$x[$i + 1] - 2/$x[$i]);
                                                 
                                        } else {
                                                 
                                                # set partial derivatives
                                                @dm = (1/$x[$i], -1/$x[$i + 1], 0, 0);
                                                 
                                        }
                                         
                                # if difference 0
                                } elsif ($dif == 0) {
                                         
                                        # if row 0
                                        if ($i == 0) {
                                                 
                                                # set partial derivatives
                                                @dm = (-3/$x[$i + 2], -2/$x[$i + 2], -1/$x[$i + 2], 0);
                                                 
                                        # if row 1
                                        } elsif ($i == 1) {
                                                 
                                                # if more than three knots
                                                if ($ix > 2) {
                                                         
                                                        # set partial derivatives
                                                        @dm = (2/$x[$i + 1] + 1/$x[$i + 2], 1/$x[$i + 1], -1/$x[$i + 2], 0);
                                                         
                                                } else {
                                                         
                                                        # set partial derivatives
                                                        @dm = (2/$x[$i + 1] + 1/$x[$i + 2], 1/$x[$i + 1], -1/$x[$i + 2], -2/$x[$i + 2] - 1/$x[$i + 1]);
                                                         
                                                }
                                                 
                                        # if next to last row
                                        } elsif ($i == $ix - 1) {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 1/$x[$i + 1], -1/$x[$i + 2], -2/$x[$i + 2] - 1/$x[$i + 1]);
                                                 
                                        # if last row
                                        } elsif ($i == $ix) {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 1/$x[$i + 1], 2/$x[$i + 1], 3/$x[$i + 1]);
                                                 
                                        } else {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 1/$x[$i + 1], -1/$x[$i + 2], 0);
                                                 
                                        }
                                         
                                # if difference 1
                                } elsif ($dif == 1) {
                                         
                                        # if row 0
                                        if ($i == 0) {
                                                 
                                                # set partial derivatives
                                                @dm = (3/$x[$i + 2] + 2/$x[$i + 3], 2/$x[$i + 2] + 1/$x[$i + 3], 1/$x[$i + 2], -1/$x[$i + 3]);
                                                 
                                        # if row 1
                                        } elsif ($i == 1) {
                                                 
                                                # if more than 3 knots
                                                if ($ix > 2) {
                                                         
                                                        # set partial derivatives
                                                        @dm = (-1/$x[$i + 2], 0, 1/$x[$i + 2], -1/$x[$i + 3]);
                                                         
                                                } else {
                                                         
                                                        # set partial derivatives
                                                        @dm = (-1/$x[$i + 2], 0, 1/$x[$i + 2], 2/$x[$i + 2]);
                                                         
                                                }
                                                 
                                        # if next to last row
                                        } elsif ($i == $ix - 1) {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 0, 1/$x[$i + 2], 2/$x[$i + 2]);
                                                 
                                        } else {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 0, 1/$x[$i + 2], -1/$x[$i + 3]);
                                                 
                                        }
                                         
                                # if difference 2
                                } elsif ($dif == 2) {
                                         
                                        # if row 0
                                        if ($i == 0) {
                                                 
                                                # set partial derivatives
                                                @dm = (-2/$x[$i + 3], -1/$x[$i + 3], 0, 1/$x[$i + 3]);
                                                 
                                        } else {
                                                 
                                                # set partial derivatives
                                                @dm = (0, 0, 0, 1/$x[$i + 3]);
                                                 
                                        }
                                 
                                }
                                 
                                # compute ∂d1/∂n
                                $dd1 = akima_dpd($m[$i], $m[$i + 1], $dm[0], $dm[1], $d1);
                                 
                                # compute ∂d2/∂n
                                $dd2 = akima_dpd($m[$i + 2], $m[$i + 3], $dm[2], $dm[3], $d2);
                                 
                                # compute ∂s/∂n
                                $derv->[$i][$j] = akima_spd($d1, $d2, $m[$i + 1], $m[$i + 2], $dd1, $dd2, $dm[1], $dm[2]);
                                 
                        }
                 
                }
                 
        }
        # set object
        $self->[5] = bless($derv, 'Math::Matrix');
}
# compute partial derivative function for _akima_para
# parameters: (d1, d2, m1, m2, ∂d1/∂n, ∂d2/∂n, ∂m1/∂n, ∂m2/∂n)
# returns; (∂s/∂n)
sub akima_spd {
HideShow 10 lines of Pod
=pod
        https://www.wolframalpha.com/input/?i=d%2Fdx+(a(x)+*+c(x)+%2B+b(x)+*+d(x))%2F(a(x)+%2B+b(x))
        d/dx((a(x) c(x) + b(x) d(x))/(a(x) + b(x))) =
        (c(x) a'(x) + a(x) c'(x) + d(x) b'(x) + b(x) d'(x))/(a(x) + b(x)) - ((a'(x) + b'(x)) (a(x) c(x) + b(x) d(x)))/(a(x) + b(x))^2
=cut
        # return partial derivative
        return(($_[2] * $_[5] + $_[1] * $_[6] + $_[3] * $_[4] + $_[0] * $_[7])/($_[1] + $_[0]) - (($_[5] + $_[4]) * ($_[1] * $_[2] + $_[0] * $_[3]))/($_[1] + $_[0])**2);
}
# compute partial derivative function for _akima_para
# parameters: (m1, m2, ∂m1/∂n, ∂m2/∂n, d)
# returns; (∂d/∂n)
sub akima_dpd {
HideShow 10 lines of Pod
=pod
        https://www.wolframalpha.com/input/?i=d%2Fdx+(sqrt((a(x)+-+b(x))%5E2+%2B+c))
        d/dx(sqrt((a(x) - b(x))^2 + c)) =
        ((a(x) - b(x)) (a'(x) - b'(x)))/sqrt((a(x) - b(x))^2 + c)
=cut
        # return partial derivative
        return(($_[0] - $_[1]) * ($_[2] - $_[3])/$_[4]);
}
# add local min/max values
# parameter: (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);
                 
                # add min/max s-values
                push(@s, map {$_ + $i} @t);
                 
                # add inner knot s-value, if derivative is 0
                push(@s, $i) if ($i && $self->[3][$i] == 0);
                 
                # 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]];
                 
        }
        # 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 ($dx, $y1, $y2, $m1, $m2);
        my ($a, $b, $c, $dscr, @t);
        # compute delta x-value (x-values : range)
        $dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
        # get endpoint values
        $y1 = $self->[2][$low];
        $y2 = $self->[2][$low + 1];
        $m1 = $self->[3][$low] * $dx;
        $m2 = $self->[3][$low + 1] * $dx;
        # 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 ($dx, $h00, $h01, $h10, $h11);
        # compute delta x-value (x-values : range)
        $dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
        # 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]/$dx + $h01 * $self->[2][$low + 1]/$dx + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1])/$dx);
}
# 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 ($dx, $tc, $ttc, $h00, $h01, $h10, $h11);
        # compute delta x-value (x-values : range)
        $dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
        # 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]/$dx + $h01 * $self->[2][$low + 1]/$dx + $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 ($dx, $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) {
                 
                # compute delta x-value (x-values : range)
                $dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
                 
                # interpolate value (Lapack XS)
                # if ICC::Support::Lapack module is loaded
                if ($lapack) {
                         
                        return(ICC::Support::Lapack::hermite($t, $self->[2][$low], $self->[2][$low + 1], $dx * $self->[3][$low], $dx * $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 * $dx * $self->[3][$low] + $h11 * $dx * $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 ($dx, $y1, $y2, $m1, $m2);
        my ($A, $B, $C, $D, $dscr, @t);
        my ($d0, $d1, $cs, $cc, @r, $ccr, $sol, $lim0, $lim1);
        # compute delta x-value (x-values : range)
        $dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
        # get endpoint values
        $y1 = $self->[2][$low];
        $y2 = $self->[2][$low + 1];
        $m1 = $self->[3][$low] * $dx;
        $m2 = $self->[3][$low + 1] * $dx;
        # 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 (abs($dscr) < 5E-15) {
                         
                        # 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 (abs($dscr) < 5E-15) {
                         
                        # if ∆0 is zero
                        if (abs($d0) < 5E-15) {
                                 
                                # 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 {$a <=> $b} grep {$_ > $lim0 && $_ < $lim1} @t);
}
# test vector values
# returns true if values are unique
# parameter: (vector)
# returns: (flag)
sub _unique {
        # get parameter
        my $vec = shift();
        # local variable
        my (%x);
        # make count hash
        for (@{$vec}) {$x{$_}++}
        # return
        return(@{$vec} == keys(%x));
}
# sort object input/output values
# parameter: (ref_to_object)
sub _sort_values {
        # get parameter
        my $self = shift();
        # local variable
        my (@p);
        # if input contains x-values
        if ($#{$self->[1]} > 1) {
                 
                # warn if vectors not same length
                ($#{$self->[1]} == $#{$self->[2]}) or carp('input and output vectors not same length');
                 
                # for each point
                for my $i (0 .. $#{$self->[1]}) {
                         
                        # copy x-y values
                        $p[$i] = [$self->[1][$i], $self->[2][$i]];
                         
                }
                 
                # sort points by input value
                @p = sort {$a->[0] <=> $b->[0]} @p;
                 
                # for each point
                for my $i (0 .. $#{$self->[1]}) {
                         
                        # copy x-y values
                        $self->[1][$i] = $p[$i]->[0];
                        $self->[2][$i] = $p[$i]->[1];
                         
                }
                 
        } else {
                 
                # if x-range is decreasing
                if ($self->[1][0] > $self->[1][-1]) {
                         
                        # reverse x and y vectors
                        @{$self->[1]} = reverse(@{$self->[1]});
                        @{$self->[2]} = reverse(@{$self->[2]});
                         
                }
                 
        }
         
}
# fix highlight region data
# corrects the print-specific
# problem of highlight dot retention
# samples below the limit (0 - 1) are fixed
# parameter: (ref_to_object, limit)
sub _fix_hl {
        # get parameter
        my ($self, $limit) = @_;
        # local variable
        my ($span, $lo, $hi, $n, @ix, @xr, @x, @y, $bern, @new);
        # compute 'output' range
        $span = ($self->[2][-1] - $self->[2][0]);
        # set 'output' limits
        $lo = $self->[2][0] + 0.002 * $span;
        $hi = $self->[2][0] + 2 * $limit * $span;
        # get upper index of 'output' array
        $n = $#{$self->[2]};
        # find samples within limits
        @ix = grep {$self->[2][$_] >= $lo && $self->[2][$_] <= $hi} (1 .. $n);
        # use at least 4 samples
        @ix = ($ix[0] .. $ix[0] + 3) if (@ix < 4);
        # if input is a range
        if (@{$self->[1]} == 2) {
                 
                # interpolate x-values
                @xr = map {(1 - $_/$n) * $self->[1][0] + ($_/$n) * $self->[1][1]} (0 .. $n);
                 
                # get 'fit' x-values
                @x = @xr[@ix];
                 
        } else {
                 
                # get 'fit' x-values
                @x = @{$self->[1]}[@ix];
                 
        }
        # get 'fit' y-values
        @y = @{$self->[2]}[@ix];
        # make 'bern' object of degree 3, fitted to x-y values (requires Lapack module)
        $bern = ICC::Support::bern->new({'input' => [$x[0], $x[-1]], 'output' => [(undef) x 4], 'fit' => [\@x, \@y]});
        # set 'output' limit
        $hi = $self->[2][0] + $limit * $span;
        # for each sample
        for my $i (reverse(0 .. $n)) {
                 
                # if 'output' value ≤ limit
                if (@new || $self->[2][$i] <= $hi) {
                         
                        # compute new 'output' value using 'bern' curve
                        push(@new, $bern->transform(@{$self->[1]} > 2 ? $self->[1][$i] : $xr[$i]));
                         
                }
                 
        }
        # splice new values, saving originals in hash
        $self->[0]{'hl_original'} = [splice(@{$self->[2]}, 0, @new, reverse(@new))];
        # save x-axis intercept (highlight dot retention value)
        $self->[0]{'hl_retention'} = $bern->inverse(0);
        # compute knot derivatives
        _objderv($self);
        # add segment min/max y-values
        _minmax($self);
}
# fix shadow region data
# corrects the print-specific problem of
# non-monotonic data is the shadow region
# samples above the limit (0 - 1) are fixed
# parameter: (ref_to_object, limit)
sub _fix_sh {
        # get parameter
        my ($self, $limit) = @_;
        # local variable
        my ($span, $lo, $n, @ix, @xr, @x, @y, $bern, @new);
        # compute 'output' range
        $span = ($self->[2][-1] - $self->[2][0]);
        # set 'output' limit
        $lo = $self->[2][-1] - 2 * (1 - $limit) * $span;
        # get upper index of 'output' array
        $n = $#{$self->[2]};
        # find samples with 'output' > limit
        @ix = grep {$self->[2][$_] >= $lo} (0 .. $n);
        # use at least 4 samples
        @ix = (($ix[-1] - 3) .. $ix[-1]) if (@ix < 4);
        # if input is a range
        if (@{$self->[1]} == 2) {
                 
                # interpolate x-values
                @xr = map {(1 - $_/$n) * $self->[1][0] + ($_/$n) * $self->[1][1]} (0 .. $n);
                 
                # get 'fit' x-values
                @x = @xr[@ix];
                 
        } else {
                 
                # get 'fit' x-values
                @x = @{$self->[1]}[@ix];
                 
        }
        # get 'fit' y-values
        @y = @{$self->[2]}[@ix];
        # make 'bern' object of degree 3, fitted to x-y values (requires Lapack module)
        $bern = ICC::Support::bern->new({'input' => [$x[0], $x[-1]], 'output' => [(undef) x 4], 'fit' => [\@x, \@y]});
        # set 'output' limit
        $lo = $self->[2][0] + $limit * $span;
        # for each sample
        for my $i (0 .. $n) {
                 
                # if 'output' value ≥ limit
                if (@new || $self->[2][$i] >= $lo) {
                         
                        # compute new 'output' value using 'bern' curve
                        push(@new, $bern->transform(@{$self->[1]} > 2 ? $self->[1][$i] : $xr[$i]));
                         
                }
                 
        }
        # splice new values, saving originals in hash
        $self->[0]{'sh_original'} = [splice(@{$self->[2]}, -@new, @new, @new)];
        # compute knot derivatives
        _objderv($self);
        # add segment min/max y-values
        _minmax($self);
}
# set object contents from parameter hash
# parameters: (ref_to_object, ref_to_parameter_hash)
sub _new_from_hash {
        # get parameters
        my ($self, $hash) = @_;
        # local variables
        my ($in, $out, $type, $damp, $fix, $limit, $pok, $fok);
        # if 'input' is defined
        if (defined($in = $hash->{'input'})) {
                 
                # if 'input' is a vector (2 or more values)
                if (ICC::Shared::is_num_vector($in) && 2 <= @{$in}) {
                         
                        # set object 'input' array
                        $self->[1] = [@{$in}];
                         
                # if 'input' is an integer > 0 (2 or more values)
                } elsif (Scalar::Util::looks_like_number($in) && $in == int($in) && $in > 0) {
                         
                        # make identity curve
                        $self->[1] = [map {$_/$in} (0 .. $in)];
                         
                } else {
                         
                        # error
                        croak("invalid 'input' values");
                         
                }
                 
        } else {
                 
                # use default
                $self->[1] = [0, 1];
                 
        }
        # if 'output' is defined
        if (defined($out = $hash->{'output'})) {
                 
                # if 'output' is a vector (2 or more values)
                if (ICC::Shared::is_num_vector($out) && 2 <= @{$out}) {
                         
                        # set object 'output' array
                        $self->[2] = [@{$out}];
                         
                # if 'output' is an integer > 0 (2 or more values)
                } elsif (Scalar::Util::looks_like_number($out) && $out == int($out) && $out > 0) {
                         
                        # make identity curve
                        $self->[2] = [map {$_/$out} (0 .. $out)];
                         
                } else {
                         
                        # error
                        croak("invalid 'output' values");
                         
                }
                 
        } else {
                 
                # use default
                $self->[2] = [0, 1];
                 
        }
        # if 'type' is defined
        if (defined($type = $hash->{'type'})) {
                 
                # verify supported type
                ($type =~ m/^(natural|akima)$/) or croak('unsupported spline type');
                 
                # set object type in hash
                $self->[0]{'type'} = $type;
                 
        }
        # if 'damp' is defined
        if (defined($damp = $hash->{'damp'})) {
                 
                # verify supported type
                (Scalar::Util::looks_like_number($damp) && $damp >= 0) or croak('invalid akima damping value');
                 
                # set object type in hash
                $self->[0]{'damp'} = $damp;
                 
        }
        # verify number of input values
        (@{$self->[1]} == 2 || @{$self->[1]} == @{$self->[2]}) or croak('wrong number of input values');
        # sort input/output values
        _sort_values($self);
        # compute knot derivatives
        _objderv($self);
        # add segment min/max y-values
        _minmax($self);
        # if 'fix_hl' is defined
        if (defined($fix = $hash->{'fix_hl'})) {
                 
                # verify polarity
                ($pok = $self->[2][-1] > $self->[2][0]) or carp("'fix_hl' requires increasing output values");
                 
                # verify fix limit (-1 to 1)
                ($fok = (Scalar::Util::looks_like_number($fix) && $fix >= -1 && $fix <= 1)) or carp("invalid 'fix_hl' limit");
                 
                # compute 'output' limit
                $limit = $self->[2][0] * (1 - $fix) + $self->[2][-1] * $fix;
                 
                # fix the data if polarity and limit are OK, and spline is non-monotonic in highlights
                _fix_hl($self, abs($fix)) if ($pok && $fok && ($fix < 0 || grep {$_ <= $limit} monotonic($self, 'output')));
                 
        }
        # if 'fix_sh' is defined
        if (defined($fix = $hash->{'fix_sh'})) {
                 
                # verify polarity
                ($pok = $self->[2][-1] > $self->[2][0]) or carp("'fix_sh' requires increasing output values");
                 
                # verify fix limit (-1 to 1)
                ($fok = (Scalar::Util::looks_like_number($fix) && $fix >= -1 && $fix <= 1)) or carp("invalid 'fix_sh' limit");
                 
                # compute 'output' limit
                $limit = $self->[2][0] * (1 - $fix) + $self->[2][-1] * $fix;
                 
                # fix the data if polarity and limit are OK, and spline is non-monotonic in shadows
                _fix_sh($self, abs($fix)) if ($pok && $fok && ($fix < 0 || grep {$_ >= $limit} monotonic($self, 'output')));
                 
        }
         
}
# 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) or 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;