diff --git a/Basic/Makefile.PL b/Basic/Makefile.PL
index f7cede635..3ed79683e 100644
--- a/Basic/Makefile.PL
+++ b/Basic/Makefile.PL
@@ -12,6 +12,8 @@ my %pm = map { my $h = '$(INST_LIBDIR)/';
 $pm{'PDLdb.pl'} = '$(INST_LIB)/PDLdb.pl';
 $pm{'threads.pm'} = '$(INST_LIB)/PDL/Parallel/threads.pm';
 $pm{'SIMD.pm'} = '$(INST_LIB)/PDL/Parallel/threads/SIMD.pm';
+$pm{'TestPDL.pm'} = '$(INST_LIB)/Test/PDL.pm';
+$pm{'TestDeepPDL.pm'} = '$(INST_LIB)/Test/Deep/PDL.pm';
 
 my %man3pods = map { my $h = '$(INST_MAN3DIR)/';
 		     $h .= 'PDL::' if $_ !~ /PDL.pm$/;
diff --git a/Basic/TestDeepPDL.pm b/Basic/TestDeepPDL.pm
new file mode 100644
index 000000000..aa300ffd4
--- /dev/null
+++ b/Basic/TestDeepPDL.pm
@@ -0,0 +1,70 @@
+use strict;
+use warnings;
+
+package Test::Deep::PDL;
+
+# ABSTRACT: Test ndarrays inside data structures with Test::Deep
+
+=for Pod::Coverage init descend diag_message renderExp renderGot
+
+=cut
+
+use Test::Deep::Cmp;
+require Test::PDL;
+
+=head1 DESCRIPTION
+
+This is just an implementation class. Look at the documentation for test_pdl()
+in L<Test::PDL>.
+
+=cut
+
+sub init
+{
+	my $self = shift;
+	my $expected = shift;
+	die "Supplied value is not an ndarray" unless eval { $expected->isa('PDL') };
+	$self->{expected} = $expected;
+}
+
+sub descend
+{
+	my $self = shift;
+	my $got = shift;
+	( my $ok, $self->data->{diag} ) = Test::PDL::eq_pdl( $got, $self->{expected} );
+	return $ok;
+}
+
+sub diag_message
+{
+	my $self = shift;
+	my $where = shift;
+	return "Comparing $where as an ndarray:\n" . $self->data->{diag};
+}
+
+sub renderExp
+{
+	my $self = shift;
+	return $self->renderGot( $self->{expected} );
+}
+
+sub renderGot
+{
+	my $self = shift;
+	my $val = shift;
+	my $fmt = '%-8T %-12D (%-5S) ';
+	return eval { $val->isa('PDL') } ? ($val->info($fmt) . $val) : ("(" . Test::Deep::render_val($val) . ")");
+}
+
+=head1 BUGS
+
+The implementation of this class depends on undocumented subroutines in
+L<Test::Deep>. This may break if L<Test::Deep> gets refactored.
+
+=head1 SEE ALSO
+
+L<PDL>, L<Test::PDL>, L<Test::Deep>
+
+=cut
+
+1;
diff --git a/Basic/TestPDL.pm b/Basic/TestPDL.pm
new file mode 100644
index 000000000..a560b699b
--- /dev/null
+++ b/Basic/TestPDL.pm
@@ -0,0 +1,449 @@
+package Test::PDL;
+
+# ABSTRACT: Test Perl Data Language arrays (a.k.a. ndarrays) for equality
+
+=head1 SYNOPSIS
+
+	use PDL;
+	use Test::More tests => 3;
+	use Test::PDL qw( is_pdl :deep );
+
+	# an example of a test that succeeds
+	$got      = sequence 5;
+	$expected = pdl( 0,1,2,3,4 );
+	is_pdl( $got, $expected, 'sequence() works as expected' );
+	#   OUTPUT:
+	# ok 1 - sequence() works as expected
+
+	# if a test fails, detailed diagnostics are printed; the output is
+	# similar to that of is() from L<Test::More>
+	$got      = pdl( 0,-1,-2,3,4 );
+	$expected = sequence 5;
+	is_pdl( $got, $expected, 'demonstrate the output of a failing test' );
+	#   OUTPUT:
+	# not ok 2 - demonstrate the output of a failing test
+	#
+	#   Failed test 'demonstrate the output of a failing test'
+	#   at maint/pod.t line 16.
+	#     2/5 values do not match
+	#          got: Double   D [5]        (P    ) [0 -1 -2 3 4]
+	#     expected: Double   D [5]        (P    ) [0 1 2 3 4]
+	# First <=5 values differ at:
+	# [
+	#  [1]
+	#  [2]
+	# ]
+	# Those 'got' values: [-1 -2]
+	# Those 'expected' values: [1 2]
+
+	# ndarrays within other data structures can be tested with Test::Deep
+	use Test::Deep qw( cmp_deeply );
+	$got      = { name => 'Histogram', data => long( 17,0,1 ) };
+	$expected = { name => 'Histogram', data => test_long( 17,0,0,1 ) };
+	cmp_deeply( $got, $expected, 'demonstrate the output of a failing deep comparison' );
+	#   OUTPUT:
+	# not ok 3 - demonstrate the output of a failing deep comparison
+	#
+	#   Failed test 'demonstrate the output of a failing deep comparison'
+	#   at maint/pod.t line 30.
+	# Comparing $data->{"data"} as an ndarray:
+	# dimensions do not match in extent
+	#    got : Long     D [3]        (P    ) [17 0 1]
+	# expect : Long     D [4]        (P    ) [17 0 0 1]
+
+=cut
+
+use strict;
+use warnings;
+use PDL::Lite;
+use PDL::Types ();
+
+use base qw( Exporter );
+our @EXPORT = qw( is_pdl );
+our @EXPORT_OK = qw( eq_pdl is_pdl test_pdl );
+our %EXPORT_TAGS = ( deep => [ qw( test_pdl ) ] );
+our $VERSION = '0.22';
+
+=head1 DESCRIPTION
+
+With Test::PDL, you can compare two ndarrays for equality. The comparison is
+performed as thoroughly as possible, comparing types, dimensions, bad value
+patterns, and finally the values themselves. The exact behaviour can be
+configured by setting certain package-wide defaults (see %DEFAULTS below), or
+by supplying options in a function call.
+Test::PDL is mostly useful in test scripts.
+
+Test::PDL is to be used with the Perl Data Language (L<PDL>).
+
+By default, Test::PDL exports only one function: is_pdl(). The other functions
+are exported on demand only. The export tag C<:deep> exports test_pdl() and one
+function for each PDL type constructor (like short(), double(), etc.), prefixed
+with C<test_>: test_short(), test_double(), ...
+
+=head1 VARIABLES
+
+=head2 %DEFAULTS
+
+The default comparison criteria used by Test::PDL can be configured by setting
+the values in the %DEFAULTS hash. This can be done directly, by addressing
+%Test::PDL::DEFAULTS directly.
+
+=over 4
+
+=item atol
+
+The absolute tolerance used to compare values. Initially set to 1e-6.
+
+=item require_equal_types
+
+If true, only ndarrays with equal type can be considered equal. If false, the
+types of the ndarrays being compared is not taken into consideration. Defaults
+to true: types must match for the comparison to succeed. If you want to
+write tests like
+
+	is_pdl( $got, pdl([ 1, 3, 5, 6 ]) );
+
+without having to worry about the type of the ndarray being exactly I<double>
+(which is the default type of the pdl() constructor), set I<require_equal_types> equal to
+0.
+
+=item rtol
+
+The relative tolerance used to compare values. Initially set to 1e-6.
+
+=back
+
+=cut
+
+our %DEFAULTS = (
+	atol                => 1e-6,
+	require_equal_types => 1,
+	rtol                => 1e-6,
+);
+
+=head1 FUNCTIONS
+
+=head2 import
+
+Custom importer that recognizes configuration defaults specified at use time, as
+in
+
+	use Test::PDL -require_equal_types => 0;
+
+=cut
+
+sub import
+{
+	my $i = 0;
+	while( $i < @_ ) {
+		if( $_[ $i ] =~ /^-/ ) {
+			my( $key, $val ) = splice @_, $i, 2;
+			$key =~ s/^-(.*)/$1/;
+			PDL::barf( "invalid name $key" ) unless grep { $key eq $_ } keys %DEFAULTS;
+			PDL::barf( "undefined value for $key" ) unless defined $val;
+			$DEFAULTS{ $key } = $val;
+		}
+		else { $i++ }
+	}
+	__PACKAGE__->export_to_level( 1, @_ );
+}
+
+=head2 is_pdl
+
+=for ref # PDL
+
+Run a test comparing an ndarray to an expected ndarray, and fail with detailed
+diagnostics if they don't compare equal.
+
+=for usage # PDL
+
+	is_pdl( $got, $expected );
+	is_pdl( $got, $expected, $test_name );
+	is_pdl( $got, $expected, { test_name => $test_name } );
+	is_pdl( $got, $expected, { atol => $absolute_tolerance, ... } );
+
+Yields ok if the first two arguments are ndarrays that compare equal, not ok if
+the ndarrays are different, or if at least one is not an ndarray. Prints a
+diagnostic when the comparison fails, with the reason and a brief printout of
+both arguments. See the documentation of eq_pdl() for the comparison
+criteria. $test_name is optional.
+
+Named after is() from L<Test::More>.
+
+=cut
+
+sub is_pdl {
+  require Test::Builder;
+  my ( $got, $expected, $arg ) = @_;
+  my $tb = Test::Builder->new;
+  $tb->croak('error in arguments: third argument is an ndarray')
+    if eval { $arg->isa('PDL') };
+  my $opt = { %DEFAULTS };
+  my $name;
+  if ($arg) {
+    if (ref $arg eq 'HASH') { $opt = { %$opt, %$arg } }
+    else { $name = $arg }
+  }
+  $name ||= $opt->{test_name} || "ndarrays are equal";
+  my ($ok, $reason, $mask) = eq_pdl($got, $expected, $opt);
+  return $tb->ok(1, $name) if $ok;
+  my $rc = $tb->ok( 0, $name );
+  my $fmt = '%-8T %-12D (%-5S) ';
+  my @mismatch;
+  if (defined $mask) {
+    my $coords = defined $mask ? $mask->not->whichND : undef;
+    $coords = $coords->slice(',0:4') if defined $coords and $coords->dim(1) > 5;
+    push @mismatch, (
+      "\nFirst <=5 values differ at:", $coords,
+      "Those 'got' values: ", $got->indexND($coords),
+      "\nThose 'expected' values: ", $expected->indexND($coords),
+    );
+  }
+  $tb->diag(
+    "    $reason\n",
+    "         got: ", eval { $got->isa('PDL')      && !$got->isnull      } ? $got->info( $fmt )      : '', $got, "\n",
+    "    expected: ", eval { $expected->isa('PDL') && !$expected->isnull } ? $expected->info( $fmt ) : '', $expected,
+    @mismatch,
+  );
+  return $rc;
+}
+
+=head2 eq_pdl
+
+=for ref # PDL
+
+Return true if two ndarrays compare equal, false otherwise. In list context,
+additionally returns a diagnostic string.
+
+=for usage # PDL
+
+	my $equal = eq_pdl( $got, $expected );
+	my $equal = eq_pdl( $got, $expected, { atol => $absolute_tolerance, ... } );
+	my( $equal, $diag ) = eq_pdl( $got, $expected );
+	my( $equal, $diag ) = eq_pdl( $got, $expected, { atol => $absolute_tolerance, ... } );
+
+eq_pdl() contains just the comparison part of is_pdl(), without the
+infrastructure required to write tests with L<Test::More>. It could be used as
+part of a larger test in which the equality of two ndarrays must be verified. By
+itself, eq_pdl() does not generate any output, so it should be safe to use
+outside test suites.
+
+In list context, eq_pdl() returns a list with three elements, the first one being
+a boolean whether the ndarrays compared equal, the second being a diagnostic
+string explaining why the comparison failed (or the empty string, if it didn't
+fail). The third is either the mask of not-equal if the values didn't
+match, or C<undef>.
+This is useful in combination with L<Test::Deep>, but might also be
+useful on its own.
+
+eq_pd() does not need L<Test::Builder>, so you can use it as part of something
+else, without side effects (like generating output).
+
+The criteria for equality are the following:
+
+=over 4
+
+=item *
+
+Both arguments must be ndarrays for the comparison to succeed. Currently, there
+is no implicit conversion from scalar to ndarray.
+
+=item *
+
+The type of both ndarrays must be equal if (and only if) I<require_equal_types> is true.
+
+=item *
+
+The number of dimensions must be equal. That is, a two-dimensional ndarray only
+compares equal with another two-dimensional ndarray.
+
+=item *
+
+The extent of the dimensions are compared one by one and must match. That is, a
+ndarray with dimensions (5,4) cannot compare equal with an ndarray of dimensions
+(5,3). Note that degenerate dimensions are not treated specially, and thus a
+ndarray with dimensions (5,4,1) is considered different from an ndarray with
+dimensions (5,4).
+
+=item *
+
+For ndarrays that conform in type and shape, the bad value pattern is examined.
+If the two ndarrays have bad values in different positions, the ndarrays are
+considered different. Note that two ndarrays may compare equal even though their
+bad flag is different, if there are no bad values.
+
+=item *
+
+And last but not least, the values themselves are examined one by one.
+As of 0.21, both integer and floating-point types are compared approximately.
+The approximate comparison is
+implemented using a combination of relative and absolute tolerances, which can
+be set by supplying an argument to C<use Test::PDL>, or by supplying an
+optional hash to this function. By default, the absolute and relative
+tolerances are both equal to 1e-6. The user can specify a pure relative
+tolerance by specifying C<atol =E<gt> 0>, and a pure absolute tolerance by
+specifying C<rtol =E<gt> 0>. If both tolerances are specified, values compare
+equal if I<either> their difference is lower than or equal to the absolute
+tolerance I<or> their relative difference (with respect to the expected
+value) is lower than or equal to the relative tolerance. For expected
+values equal to zero, relative differences (with respect to the expected
+value) make no sense, and the use of combined absolute and relative
+tolerances is recommended.
+
+=back
+
+=cut
+
+sub eq_pdl {
+  my ($got, $expected, $arg) = @_;
+  my $opt = { %DEFAULTS, ref $arg eq 'HASH' ? %$arg : () };
+  PDL::barf( 'need an absolute or a relative tolerance, or both' ) unless defined $opt->{atol} || defined $opt->{rtol};
+  $opt->{atol} //= 0;
+  $opt->{rtol} //= 0;
+  PDL::barf('absolute tolerance cannot be negative') if $opt->{atol} < 0;
+  PDL::barf('relative tolerance cannot be negative') if $opt->{rtol} < 0;
+  return wantarray ? (0, 'received value is not an ndarray', undef) : 0
+    if !eval { $got->isa('PDL') };
+  return wantarray ? (0, 'expected value is not an ndarray', undef) : 0
+    if !eval { $expected->isa('PDL') };
+  return wantarray ? (0, 'types do not match (\'require_equal_types\' is true)', undef) : 0
+    if $opt->{require_equal_types} && $got->type != $expected->type;
+  my @got_dims = $got->dims;
+  my @exp_dims = $expected->dims;
+  return wantarray ? (0, 'dimensions do not match in number', undef) : 0
+    if @got_dims != @exp_dims;
+  while (@got_dims) {
+    return wantarray ? (0, 'dimensions do not match in extent', undef) : 0
+      if shift(@got_dims) != shift(@exp_dims);
+  }
+  return wantarray ? (1, '', undef) : 1
+    if $got->isempty and $expected->isempty;
+  # both are now non-empty
+  my $res = PDL::Primitive::approx_artol( $got, $expected, @$opt{qw(atol rtol)} );
+  return wantarray ? (1, '', undef) : 1 if $res->all;
+  my $exp_nelem = $expected->nelem;
+  my $reason = ($exp_nelem-$res->sum)."/$exp_nelem values do not match";
+  return wantarray ? (0, $reason, $res) : 0;
+}
+
+=head2 test_pdl
+
+=for ref # PDL
+
+Special comparison to be used in conjunction with L<Test::Deep> to test ndarrays
+inside data structures.
+
+=for usage # PDL
+
+	my $expected = { ..., some_field => test_pdl( 1,2,-7 ), ... };
+	my $expected = [ ..., test_short( 1,2,-7 ), ... ];
+
+Suppose you want to compare data structures that happen to contain ndarrays. You
+use is_deeply() (from L<Test::More>) or cmp_deeply() (from L<Test::Deep>) to
+compare the structures element by element. Unfortunately, you cannot just write
+
+	my $got = my_sub( ... );
+	my $expected = {
+		...,
+		some_field => pdl( ... ),
+		...
+	};
+	is_deeply $got, $expected;
+
+Neither does cmp_deeply() work in the same situation. is_deeply() tries to
+compare the ndarrays using the (overloaded) C<==> comparison operator, which
+doesn't work. It simply dies with an error message saying that multidimensional
+ndarrays cannot be compared, whereas cmp_deeply() performs only a shallow
+comparison of the references.
+
+What you need is a special comparison, which is provided by this function, to
+be used with cmp_deeply(). You need to rewrite $expected as follows
+
+	my $expected = {
+		...,
+		some_field => test_pdl( ... ),
+		...
+	};
+	cmp_deeply $got, $expected;
+
+Note that you need to write test_pdl() instead of pdl(). You could achieve the
+same thing with
+
+	my $expected = {
+		...,
+		some_field => code( sub { eq_pdl( shift, pdl( ... ) ) } ),
+		...
+	};
+
+but the diagnostics provided by test_pdl() are better, and it's easier to use.
+test_pdl() accepts the same arguments as the PDL constructor pdl() does. If you
+need to compare an ndarray with a type different from the default type, use one
+of the provided test_byte(), test_short(), test_long(), etc.:
+
+	my $expected = { data => test_short( -4,-9,13 ) };
+
+If you need to manipulate the expected value, you should keep in mind that the
+return value of test_pdl() and the like are not ndarrays. Therefore, in-place
+modification of the expected value won't work:
+
+	my $expected = { data => test_short( -99,-9,13 )->inplace->setvaltobad( -99 ) }; # won't work!
+
+You should rather do
+
+	my $expected = { data => test_pdl( short(-99,-9,13)->inplace->setvaltobad(-99) ) };
+
+test_pdl() will correctly set the type of the expected value to I<short> in the
+above example.
+
+=cut
+
+sub test_pdl
+{
+	require Test::Deep::PDL;
+	my $expected = PDL::Core::pdl( @_ );
+	return Test::Deep::PDL->new( $expected );
+}
+
+=for Pod::Coverage test_anyval test_byte test_short test_ushort test_long
+test_indx test_longlong test_float test_double test_cfloat test_cdouble
+test_cldouble test_ldouble test_sbyte test_ulong test_ulonglong
+
+=cut
+
+for my $type ( PDL::Types::types ) {
+	my $sub = sub {
+		require Test::Deep::PDL;
+		my $expected = PDL::convert(
+			PDL::Core::alltopdl( 'PDL', scalar(@_) > 1 ? [@_] : shift ),
+			$type->numval
+		);
+		return Test::Deep::PDL->new( $expected );
+	};
+	my $sub_name = 'test_' . $type->convertfunc;
+	{
+		no strict 'refs';
+		*$sub_name = $sub;
+	}
+	push @EXPORT_OK, $sub_name;
+	push @{ $EXPORT_TAGS{deep} }, $sub_name;
+}
+
+=head1 BUGS
+
+None reported so far.
+
+=head1 SEE ALSO
+
+L<PDL>, L<Test::More>, L<Test::Deep>, L<Test::PDL::Deep>
+
+=head1 ACKNOWLEDGMENTS
+
+Thanks to PDL Porters Joel Berger, Chris Marshall, and David Mertens for
+feedback and improvements.
+
+Thanks to Ed J, Zakariyya Mughal, and Diab Jerius for feedback, improvements,
+maintenance of the code, and encouragement!
+
+=cut
+
+1;
diff --git a/Basic/t/tp-deep.t b/Basic/t/tp-deep.t
new file mode 100644
index 000000000..ec253aea0
--- /dev/null
+++ b/Basic/t/tp-deep.t
@@ -0,0 +1,179 @@
+use strict;
+use warnings;
+use Test::More 0.88;
+use Test::Deep qw( cmp_deeply code );
+use Test::PDL qw( :deep eq_pdl );
+use Test::Builder::Tester;
+use Test::Exception;
+use PDL;
+use PDL::Types;
+
+my @types = PDL::Types::types;
+
+isa_ok test_pdl( 1,2,3 ), 'Test::Deep::PDL';
+for my $type ( @types ) {
+	no strict 'refs';
+	my $sub = "test_$type";
+	isa_ok $sub->( 1,2,3 ), 'Test::Deep::PDL';
+}
+
+{
+	my $pdl1     = pdl( 1,2,3.13 );
+	my $got      = { name => 'Histogram', data => $pdl1 };
+	my $pdl2     = pdl( 1,2,3.13 );
+	my $expected = { name => 'Histogram', data => $pdl2 };
+	throws_ok { ok $pdl1 == $pdl2 } qr/multielement ndarray in conditional expression /, '== dies with an error message';
+	throws_ok { is $pdl1, $pdl2 } qr/multielement ndarray in conditional expression /, 'is() dies with an error message';
+}
+
+{
+	my $pdl      = pdl( 1,2,3.13 );
+	my $got      = { name => 'Histogram', data => $pdl };
+	my $expected = { name => 'Histogram', data => $pdl };
+	throws_ok { ok $pdl == $pdl } qr/^multielement ndarray in conditional expression /, 'even shallow reference comparisons do not work with ==';
+}
+
+{
+	my $pdl      = pdl( 1,2,3.13 );
+	my $got      = { name => 'Histogram', data => $pdl };
+	my $expected = { name => 'Histogram', data => $pdl };
+	test_out 'ok 1';
+	cmp_deeply $got, $expected;
+	test_test 'cmp_deeply() without test_pdl() performs only shallow reference comparison';
+}
+
+{
+	my $pdl      = pdl( 1,2,3.13 );
+	my $got      = { name => 'Histogram', data => $pdl };
+	my $expected = { name => 'Histogram', data => $pdl->copy };
+	test_out 'not ok 1';
+	test_fail +4;
+	test_diag 'Compared ${$data->{"data"}}';
+	test_err  "/#    got : '-?\\d+'/",
+		  "/# expect : '-?\\d+'/";
+	cmp_deeply $got, $expected;
+	test_test 'but shallow reference comparison is not powerful enough';
+}
+
+=pod
+
+The following test code may be a bit hard to follow. We're basically trying to
+ensure that
+
+	my @vals     = ( ... );
+	my $got      = { data => long( @vals ) };
+	my $expected = { data => test_long( @vals ) };
+	cmp_deeply $got, $expected;
+
+passes for every conceivable type of ndarray (not only I<long>), and for
+different sets of values @vals, some of which may contain bad values. We also
+test that
+
+	{ data => test_long( @vals ) }
+	{ data => test_pdl( long(@vals) ) }
+	{ data => code( sub { eq_pdl shift, long(@vals) } ) }
+
+yield the same results.
+
+=cut
+
+for my $vals ( [ 0 ], [ 2,3,0,1,99 ], [ 99,99,99 ] ) {
+	for my $type ( @types ) {
+		my @vals      = @$vals;
+		my $ctor      = do { local *slot = $PDL::{ $type }; *slot{CODE} };
+		my $tester    = do { local *slot = $Test::PDL::{ 'test_' . $type }; *slot{CODE} };
+		my $pdl       = $ctor->( @vals )->inplace->setvaltobad( 99 );
+		note 'test with pdl = ', $pdl->info('%-8T'), ' ', $pdl;
+		my $got       = { data => $pdl };
+		my $expected1 = { data => $tester->( @vals ) };
+		$expected1->{data}->{expected}->inplace->setvaltobad( 99 );
+		test_out 'ok 1';
+		cmp_deeply $got, $expected1;
+		test_test "$type succeeds when it should succeed, with ndarray supplied as values (@vals)";
+		my $expected2 = { data => test_pdl( $pdl ) };
+		test_out 'ok 1';
+		cmp_deeply $got, $expected2;
+		test_test "... ($type) also when ndarray is supplied directly (@vals)";
+		my $expected3 = { data => code( sub { eq_pdl shift, $pdl } ) };
+		test_out 'ok 1';
+		cmp_deeply $got, $expected3;
+		test_test "... ($type) and it's the same thing as using code() (@vals)";
+	}
+}
+
+{
+	my $pdl1 = 2;
+	my $pdl2 = pdl( 3,4,9.999 );
+	ok !eq_pdl( $pdl1, $pdl2 ), 'ndarrays are unequal to begin with';
+	my $got = { data => $pdl1 };
+	my $expected = { data => test_pdl( $pdl2 ) };
+	test_out 'not ok 1';
+	test_fail +5;
+	test_diag 'Comparing $data->{"data"} as an ndarray:',
+		  'received value is not an ndarray';
+	test_err  "/#    got : \\('2'\\)/",
+		  '/# expect : Double\s+D\s+\[3\].*/';
+	cmp_deeply $got, $expected;
+	test_test 'fails with correct message and diagnostics when received value is not an ndarray';
+	test_out 'not ok 1';
+	test_fail +6;
+	test_diag 'Ran coderef at $data->{"data"} on';
+	test_err  '/#?\s*/';
+	test_diag "'2'",
+		  'and it said',
+		  'received value is not an ndarray';
+	cmp_deeply $got, { data => code( sub { eq_pdl shift, $pdl2 } ) };
+	test_test '... but the diagnostics are better than with code()';
+}
+
+{
+	my $pdl1 = pdl( 3,4,9.999 );
+	my $pdl2 = pdl( 3,4,10 );
+	ok !eq_pdl( $pdl1, $pdl2 ), 'ndarrays are unequal to begin with';
+	my $got = { data => $pdl1 };
+	my $expected = { data => test_pdl( $pdl2 ) };
+	test_out 'not ok 1';
+	test_fail +5;
+	test_diag 'Comparing $data->{"data"} as an ndarray:',
+		  '1/3 values do not match';
+	test_err  '/#    got : Double\s+D\s+\[3\].*/',
+		  '/# expect : Double\s+D\s+\[3\].*/';
+	cmp_deeply $got, $expected;
+	test_test 'fails with correct message and diagnostics on value mismatch';
+	test_out 'not ok 1';
+	test_fail +6;
+	test_diag 'Ran coderef at $data->{"data"} on';
+	test_err  '/#?\s*/',
+		  '/# PDL=SCALAR\(0x[0-9A-Fa-f]+\)/';
+	test_diag 'and it said',
+		  '1/3 values do not match';
+	cmp_deeply $got, { data => code( sub { eq_pdl shift, $pdl2 } ) };
+	test_test '... but the diagnostics are better than with code()';
+}
+
+{
+	my $pdl1 = short( 3,4,-6 );
+	my $pdl2 = long( 3,4,10 );
+	ok !eq_pdl( $pdl1, $pdl2 ), 'ndarrays are unequal to begin with';
+	my $got = { data => $pdl1 };
+	my $expected = { data => test_pdl( $pdl2 ) };
+	test_out 'not ok 1';
+	test_fail +5;
+	test_diag 'Comparing $data->{"data"} as an ndarray:',
+		  'types do not match (\'require_equal_types\' is true)';
+	test_err  '/#    got : Short\s+D\s+\[3\].*/',
+		  '/# expect : Long\s+D\s+\[3\].*/';
+	cmp_deeply $got, $expected;
+	test_test 'fails with correct message and diagnostics on type mismatch';
+	test_out 'not ok 1';
+	test_fail +6;
+	test_diag 'Ran coderef at $data->{"data"} on';
+	test_err  '/#?\s*/',
+		  '/# PDL=SCALAR\(0x[0-9A-Fa-f]+\)/';
+	test_diag 'and it said',
+		  'types do not match (\'require_equal_types\' is true)';
+	cmp_deeply $got, { data => code( sub { eq_pdl shift, $pdl2 } ) };
+	test_test '... but the diagnostics are better than with code()';
+}
+
+done_testing;
diff --git a/Basic/t/tp-eq_pdl.t b/Basic/t/tp-eq_pdl.t
new file mode 100644
index 000000000..4ccf11ece
--- /dev/null
+++ b/Basic/t/tp-eq_pdl.t
@@ -0,0 +1,319 @@
+use strict;
+use warnings;
+use Test::More;
+use PDL;
+use Test::PDL qw( eq_pdl );
+my @warns; $SIG{__WARN__} = sub {push @warns, @_};
+
+sub run_eq_pdl
+{
+	my $scalar = eq_pdl(@_);
+	my @list = eq_pdl(@_);
+	cmp_ok(scalar @list, '==', 3, 'eq_pdl() returns a list with two elements in list context');
+	cmp_ok $scalar, '==', $list[0], '... and first element matches the return value in scalar context';
+	return @list;
+}
+
+my $values_not_match = qr/values do not match/;
+
+my ( $got, $expected, $ok, $diag );
+
+( $ok, $diag ) = run_eq_pdl();
+ok !$ok, 'rejects missing arguments';
+is $diag, 'received value is not an ndarray';
+
+$got = pdl( 9,-10 );
+( $ok, $diag ) = run_eq_pdl( $got );
+ok !$ok, 'rejects missing arguments';
+is $diag, 'expected value is not an ndarray';
+
+$expected = 3;
+$got = 4;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'rejects non-ndarray arguments';
+is $diag, 'received value is not an ndarray';
+
+$expected = 3;
+$got = long( 3,4 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'rejects non-ndarray arguments';
+is $diag, 'expected value is not an ndarray';
+
+$expected = short( 1,2 );
+$got = -2;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'rejects non-ndarray arguments';
+is $diag, 'received value is not an ndarray';
+
+$expected = long( 3,4 );
+$got = pdl( 3,4 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { require_equal_types => 0 } );
+ok $ok, 'all else being equal, compares equal on differing types when \'require_equal_types\' is false';
+is $diag, '';
+
+$expected = long( 3,4 );
+$got = pdl( 3,4 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { require_equal_types => 1 } );
+ok !$ok, 'catches type mismatch, but only when \'require_equal_types\' is true';
+is $diag, 'types do not match (\'require_equal_types\' is true)';
+
+$expected = long( 3 );
+$got = long( 3,4 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'catches dimensions mismatches (number of dimensions)';
+is $diag, 'dimensions do not match in number';
+
+$expected = zeroes( double, 3,4 );
+$got = zeroes( double, 3,4,1 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'does not treat degenerate dimensions specially';
+is $diag, 'dimensions do not match in number';
+
+$expected = long( [ [3,4],[1,2] ] );
+$got = long( [ [3,4,5], [1,2,3] ] );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'catches dimensions mismatches (extent of dimensions)';
+is $diag, 'dimensions do not match in extent';
+
+$expected = long( 4,5,6,-1,8,9 )->inplace->setvaltobad( -1 );
+$got = long( 4,5,6,7,-1,9 )->inplace->setvaltobad( -1 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'catches bad value pattern mismatch';
+like $diag, $values_not_match;
+
+$expected = long( 4,5,6,7,8,9 );
+$got = long( 4,5,6,7,-8,9 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'catches value mismatches for integer data';
+like $diag, $values_not_match;
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,-8,9 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'catches value mismatches for floating-point data';
+like $diag, $values_not_match;
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.001,9 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'approximate comparison for floating-point data fails correctly at documented default tolerance of 1e-6';
+like $diag, $values_not_match;
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.0000001,9 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok $ok, 'approximate comparison for floating-point data succeeds correctly at documented default tolerance of 1e-6';
+is $diag, '';
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.001,9 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-2 } );
+ok $ok, 'approximate comparison for floating-point data succeeds correctly at user-specified tolerance of 1e-2';
+is $diag, '';
+
+$expected = pdl( 0,1,2,3,4 );
+$got = sequence 5;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok $ok, 'succeeds when it should succeed';
+is $diag, '';
+
+$expected = null;
+$got = null;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok $ok, 'null == null';
+is $diag, '';
+
+$got = zeroes(0,3);
+$expected = zeroes(0,2);
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'differently-shaped empties are different';
+is $diag, 'dimensions do not match in extent';
+
+$got = zeroes(0);
+$expected = null;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok $ok, 'null == empty';
+is $diag, '';
+
+$expected = null;
+$got = pdl( 1,2,3 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'pdl( ... ) != null';
+is $diag, 'dimensions do not match in extent';
+
+$expected = pdl( 1,2,3 );
+$got = null;
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok !$ok, 'null != pdl( ... )';
+is $diag, 'dimensions do not match in extent';
+
+################################################################################
+note 'mixed-type comparisons';
+
+$expected = double( 0,1,2.001,3,4 );
+$got = long( 0,1,2,3,4 );
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+ok $ok, 'succeeds correctly for long/double';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+ok !$ok, 'fails correctly for long/double';
+like $diag, $values_not_match;
+
+$expected = short( 0,1,2,3,4 );
+$got = float( 0,1,2.001,3,4 );
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+ok $ok, 'succeeds correctly for float/short';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+ok !$ok, 'fails correctly for float/short';
+like $diag, $values_not_match;
+
+$expected = float( 0,-1,2.001,3,49999.998 );
+$got = double( 0,-0.9999,1.999,3,49999.999 );
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+ok $ok, 'succeeds correctly for double/float';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+ok !$ok, 'fails correctly for double/float';
+like $diag, $values_not_match;
+
+################################################################################
+note 'tests with values of significantly different magnitudes, no zeroes';
+$expected = double( 1e+3, 1, 1e-3 );
+$got = double( 1.001e+3, 1.001, 1.001e-3 );
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 0.999, rtol => 0 } );
+ok !$ok, 'still fails with an absolute tolerance of 0.999';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1, rtol => 0 } );
+ok $ok, 'passes with an absolute tolerance of 1';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-4, rtol => 0 } );
+ok !$ok, 'fails for case with different magnitudes and pure absolute tolerance of 1e-4';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-2, rtol => 0 } );
+ok !$ok, 'still fails with an absolute tolerance of 1e-2';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 10, rtol => 0 } );
+ok $ok, 'needs an absolute tolerance of 10 to pass';
+is $diag, '';
+
+# notice the other values that are way off ...
+( $ok, $diag ) = run_eq_pdl( double( 1.001e+3, 9, 9 ), $expected, { atol => 10, rtol => 0 } );
+ok $ok, '... and this leads to large errors in the smaller components';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 0, rtol => 1e-4 } );
+ok !$ok, 'should not pass with a pure relative tolerance of 1e-4';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 0, rtol => 1e-2 } );
+ok $ok, 'but passes with a pure relative tolerance of 1e-2';
+is $diag, '';
+
+################################################################################
+note 'tests with values of significantly different magnitudes, with zeroes';
+$expected = double( 1e+3, 1, 1e-9, 0 );
+$got = double( 1.00001e+3, .99999, 1.00001e-9, 1e-5 );
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 0, rtol => 1e-6 } );
+ok !$ok, 'fails at pure relative tolerance of 1e-6';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 0, rtol => 1e-4 } );
+ok !$ok, 'but also fails at pure relative tolerance of 1e-4';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( $got, $expected, { atol => 1e-4, rtol => 1e-4 } );
+ok $ok, 'needs both absolute and relative tolerances to pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double( 1e+3, 1, 0.001, 0 ), $expected, { atol => 1e-4, rtol => 1e-4 } );
+ok !$ok, 'combination of relative and absolute tolerances avoids large relative errors in small components'; # (provided atol is not too high)
+like $diag, $values_not_match;
+
+################################################################################
+note 'test perfect equality';
+
+( $ok, $diag ) = run_eq_pdl( pdl(1), pdl(1), { atol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( pdl(0), pdl(0), { atol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( pdl(0), pdl(0), { rtol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( pdl(-1e20), pdl(-1e20), { atol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( pdl(-1e-20), pdl(-1e-20), { atol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double( 0,0,0 ), double( 0,0,0 ), { atol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double( 0,0,0 ), double( 0,0,0 ), { rtol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double( 0,0,0 ), double( 0,0,0 ), { atol => 1e-10, rtol => 1e-10 } );
+ok $ok, 'perfectly equal ndarrays should always pass';
+is $diag, '';
+
+################################################################################
+note 'test tolerance "sharpness"';
+
+( $ok, $diag ) = run_eq_pdl( double(1.99), double(1), { atol => 0, rtol => 1 } );
+ok $ok, 'passes correctly within tolerance';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double(1.99), double(1), { atol => 1, rtol => 0 } );
+ok $ok, 'passes correctly within tolerance';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double(1.99), double(1), { atol => 1, rtol => 1 } );
+ok $ok, 'passes correctly within tolerance';
+is $diag, '';
+
+( $ok, $diag ) = run_eq_pdl( double(2.01), double(1), { atol => 0, rtol => 1 } );
+ok !$ok, 'fails correctly just outside of tolerance';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( double(2.01), double(1), { atol => 1, rtol => 0 } );
+ok !$ok, 'fails correctly just outside of tolerance';
+like $diag, $values_not_match;
+
+( $ok, $diag ) = run_eq_pdl( double(2.01), double(1), { atol => 1, rtol => 1 } );
+ok !$ok, 'combined tolerances should not yield a larger comparison margin';
+like $diag, $values_not_match;
+
+################################################################################
+note 'miscellaneous';
+
+$expected = long( 4,5,6,7,8,9 );
+$expected->badflag( 1 );
+$got = long( 4,5,6,7,8,9 );
+$got->badflag( 0 );
+( $ok, $diag ) = run_eq_pdl( $got, $expected );
+ok $ok, "isn't fooled by differing badflags";
+is $diag, '';
+
+is "@warns", "", "no warnings";
+done_testing;
diff --git a/Basic/t/tp-eq_pdl_clean.t b/Basic/t/tp-eq_pdl_clean.t
new file mode 100644
index 000000000..842c386b2
--- /dev/null
+++ b/Basic/t/tp-eq_pdl_clean.t
@@ -0,0 +1,48 @@
+use strict;
+use warnings;
+use Test::More;
+
+eval { require Capture::Tiny };
+plan skip_all => 'Capture::Tiny not found' if $@;
+
+my $expected_stderr = '';
+
+# first capture the output of a program that does *not* call eq_pdl()
+# this is to work around warnings that might be emitted by other modules (e.g.,
+# File::Map on some platforms complaining about the :all tag)
+{
+	my $rc;
+	my( $stdout, $stderr ) = Capture::Tiny::capture( sub {
+			my @cmd = ( $^X, '-Ilib', '-MTest::PDL=eq_pdl', '-e1' );
+			$rc = system @cmd;
+		} );
+	cmp_ok $rc, '==', 0, 'system() succeeded';
+	is $stdout, '', 'no output on stdout';
+	$expected_stderr = $stderr;
+}
+
+# test that eq_pdl() doesn't produce any output so it can safely be used in non-test code
+{
+	my $rc;
+	my( $stdout, $stderr ) = Capture::Tiny::capture( sub {
+			my @cmd = ( $^X, '-Ilib', '-MTest::PDL=eq_pdl', '-e', 'scalar eq_pdl(3,4)' );
+			$rc = system @cmd;
+		} );
+	cmp_ok $rc, '==', 0, 'system() succeeded';
+	is $stdout, '', 'eq_pdl() does not produce output on stdout';
+	is $stderr, $expected_stderr, 'eq_pdl() does not produce output on stderr';
+}
+
+# test that eq_pdl() doesn't produce any output in list context so it can safely be used in non-test code
+{
+	my $rc;
+	my( $stdout, $stderr ) = Capture::Tiny::capture( sub {
+			my @cmd = ( $^X, '-Ilib', '-MTest::PDL=eq_pdl', '-e', '(eq_pdl(3,4))' );
+			$rc = system @cmd;
+		} );
+	cmp_ok $rc, '==', 0, 'system() succeeded';
+	is $stdout, '', 'eq_pdl() does not produce output on stdout in list context';
+	is $stderr, $expected_stderr, 'eq_pdl() does not produce output on stderr in list context';
+}
+
+done_testing;
diff --git a/Basic/t/tp-import_options.t b/Basic/t/tp-import_options.t
new file mode 100644
index 000000000..3238da9d1
--- /dev/null
+++ b/Basic/t/tp-import_options.t
@@ -0,0 +1,92 @@
+use strict;
+use warnings;
+use Test::More;
+use Test::Deep;
+use Test::Exception;
+use PDL::Types 'types';
+my @warns; $SIG{__WARN__} = sub {push @warns, @_};
+
+is "@warns", "", "no warnings";
+require Test::PDL;
+@warns = ();
+
+# we should start out without an 'is_pdl' function
+ok ! __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL '';
+package t1;
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 1,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok ! __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL;
+package t2;
+Test::PDL->import();
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 1,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL -require_equal_types => 0;
+package t3;
+Test::PDL->import( -require_equal_types => 0 );
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 0,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+$Test::PDL::DEFAULTS{require_equal_types} = 1; # explicitly reset so no need reload
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL -atol => 1e-8;
+package t4;
+Test::PDL->import( -atol => 1e-8 );
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-8 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 1,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL -atol => 1e-8, -require_equal_types => 0, 'is_pdl';
+package t5;
+Test::PDL->import( -atol => 1e-8, -require_equal_types => 0, 'is_pdl' );
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-8 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 0,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# "reset"
+package t5;
+Test::PDL->import( -atol => 1e-6, -require_equal_types => 1, -rtol => 1e-6 );
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 1,
+	rtol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL -rtol => 1e-8;
+package t5;
+Test::PDL->import( -rtol => 1e-8 );
+::cmp_deeply \%Test::PDL::DEFAULTS, {
+	atol                => ::code( sub { abs( $_[0]/1e-6 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+	require_equal_types => 1,
+	rtol                => ::code( sub { abs( $_[0]/1e-8 - 1 ) < 1e-6 ? 1 : ( 0, 'tolerance beyond specified value' ) } ),
+};
+::ok __PACKAGE__->can( 'is_pdl' );
+
+# use Test::PDL -whatever => 42;
+package t6;
+::throws_ok { Test::PDL->import( -whatEver => 42 ) } qr/\binvalid name whatEver\b/;
+::ok ! __PACKAGE__->can( 'is_pdl' );
+
+::is "@warns", "", "no warnings";
+::done_testing;
diff --git a/Basic/t/tp-is_pdl.t b/Basic/t/tp-is_pdl.t
new file mode 100644
index 000000000..097910f9a
--- /dev/null
+++ b/Basic/t/tp-is_pdl.t
@@ -0,0 +1,216 @@
+use strict;
+use warnings;
+
+use Test::More;
+use Test::Builder::Tester;
+use Test::Exception;
+use PDL;
+use Test::PDL;
+my @warns; $SIG{__WARN__} = sub {push @warns, @_};
+
+my ( $got, $expected );
+
+my $values_not_match = '/#\s+\d+\/\d+\s+values do not match\n(.|\n)*/';
+
+$expected = 3;
+$got = long( 3,4 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+expected value is not an ndarray\n(.|\n)*/' );
+is_pdl( $got, $expected );
+test_test( 'rejects non-ndarray arguments' );
+
+$expected = short( 1,2 );
+$got = -2;
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+received value is not an ndarray\n(.|\n)*/' );
+is_pdl( $got, $expected );
+test_test( 'rejects non-ndarray arguments' );
+
+$expected = long( 3,4 );
+$got = pdl( 3,4 );
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected, { require_equal_types => 0 } );
+test_test( 'all else being equal, compares equal on differing types when \'require_equal_types\' is false' );
+
+$expected = long( 3,4 );
+$got = pdl( 3,4 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+types do not match \([\']require_equal_types[\'] is true\)\n(.|\n)*/' );
+is_pdl( $got, $expected, { require_equal_types => 1 } );
+test_test( 'catches type mismatch, but only when \'require_equal_types\' is true' );
+
+$expected = long( 3 );
+$got = long( 3,4 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+dimensions do not match in number\n(.|\n)*/' );
+is_pdl( $got, $expected );
+test_test( 'catches dimensions mismatches (number of dimensions)' );
+
+$expected = zeroes( double, 3,4 );
+$got = zeroes( double, 3,4,1 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+dimensions do not match in number\n(.|\n)*/' );
+is_pdl( $got, $expected );
+test_test( 'does not treat degenerate dimensions specially' );
+
+$expected = long( [ [3,4],[1,2] ] );
+$got = long( [ [3,4,5], [1,2,3] ] );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( '/#\s+dimensions do not match in extent\n(.|\n)*/' );
+is_pdl( $got, $expected );
+test_test( 'catches dimensions mismatches (extent of dimensions)' );
+
+$expected = long( 4,5,6,-1,8,9 )->inplace->setvaltobad( -1 );
+$got = long( 4,5,6,7,-1,9 )->inplace->setvaltobad( -1 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected );
+test_test( 'catches bad value pattern mismatch' );
+
+$expected = long( 4,5,6,7,8,9 );
+$got = long( 4,5,6,7,-8,9 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected );
+test_test( 'catches value mismatches for integer data' );
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,-8,9 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected );
+test_test( 'catches value mismatches for floating-point data' );
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.001,9 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected );
+test_test( 'approximate comparison for floating-point data fails correctly at documented default tolerance of 1e-6' );
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.0000001,9 );
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected );
+test_test( 'approximate comparison for floating-point data succeeds correctly at documented default tolerance of 1e-6' );
+
+$expected = pdl( 4,5,6,7,8,9 );
+$got = pdl( 4,5,6,7,8.001,9 );
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected, { atol => 1e-2 } );
+test_test( 'approximate comparison for floating-point data succeeds correctly at user-specified tolerance of 1e-2' );
+
+$expected = pdl( 0,1,2,3,4 );
+$got = sequence 5;
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected );
+test_test( 'succeeds when it should succeed' );
+
+$expected = null;
+$got = null;
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected );
+test_test( 'null == null' );
+
+$expected = null;
+$got = pdl( 1,2,3 );
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +4 );
+test_err( '/#\s+dimensions do not match in extent/' );
+test_err( '/#\s+got:.*/' );
+test_err( '/#\s+expected:\s+Null/' );
+is_pdl( $got, $expected, { require_equal_types => 0 } );
+test_test( 'pdl( ... ) != null' );
+
+$expected = pdl( 1,2,3 );
+$got = null;
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +4 );
+test_err( '/#\s+dimensions do not match in extent/' );
+test_err( '/#\s+got:\s+Null/' );
+test_err( '/#\s+expected:.*/' );
+is_pdl( $got, $expected, { require_equal_types => 0 } );
+test_test( 'null != pdl( ... )' );
+
+note 'mixed-type comparisons';
+
+$expected = double( 0,1,2.001,3,4 );
+$got = long( 0,1,2,3,4 );
+
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+test_test( 'succeeds correctly for long/double' );
+
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+test_test( 'fails correctly for long/double' );
+
+$expected = short( 0,1,2,3,4 );
+$got = float( 0,1,2.001,3,4 );
+
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+test_test( 'succeeds correctly for float/short' );
+
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+test_test( 'fails correctly for float/short' );
+
+$expected = float( 0,-1,2.001,3,49999.998 );
+$got = double( 0,-0.9999,1.999,3,49999.999 );
+
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected, { atol => 1e-2, require_equal_types => 0 } );
+test_test( 'succeeds correctly for double/float' );
+
+test_out( "not ok 1 - ndarrays are equal" );
+test_fail( +2 );
+test_err( $values_not_match );
+is_pdl( $got, $expected, { atol => 1e-6, require_equal_types => 0 } );
+test_test( 'fails correctly for double/float' );
+
+note 'miscellaneous';
+
+$expected = pdl( 0,1,2,3,4 );
+$got = sequence 5;
+test_out( "ok 1 - insert custom test name" );
+is_pdl( $got, $expected, 'insert custom test name' );
+test_test( 'custom test name is displayed correctly' );
+
+test_out( "ok 1 - insert custom test name" );
+is_pdl( $got, $expected, { test_name => 'insert custom test name' } );
+test_test( 'custom test name is also displayed correctly when supplied as an option hash' );
+
+# Although the next test may appear strange, the case it tests can be produced
+# by the following test:
+#	is_pdl hist( pdl(2,3,4,5) ), pdl(1,1,1,1);
+# Since hist() returns two ndarrays in list context, the expected ndarray ends up
+# as the third argument. Since this is probably not what the user intended, an
+# error is raised.
+throws_ok { is_pdl( $got, $expected, pdl(1,1,1,1) ) }
+	qr/^error in arguments: third argument is an ndarray at /, 'third argument is an ndarray';
+
+$expected = long( 4,5,6,7,8,9 );
+$expected->badflag( 1 );
+$got = long( 4,5,6,7,8,9 );
+$got->badflag( 0 );
+test_out( "ok 1 - ndarrays are equal" );
+is_pdl( $got, $expected );
+test_test( "isn't fooled by differing badflags" );
+
+is "@warns", "", "no warnings";
+done_testing;
diff --git a/Changes b/Changes
index 5389b8437..ad2c8fd9b 100644
--- a/Changes
+++ b/Changes
@@ -21,6 +21,7 @@
 - removed PDL::Config mechanism, leave vestigial file for back-compat and use env vars for libs/incs
 - incorporate PDL::Parallel::threads
 - Ufunc::firstnonzeroover added - thanks @guillepo for inspiration
+- incorporate Test::PDL
 
 2.093 2024-09-29
 - PDL.set_datatype now doesn't physicalise input, PDL.convert_type does
diff --git a/MANIFEST b/MANIFEST
index e6489a752..da6a8f592 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -272,7 +272,14 @@ Basic/t/slice.t
 Basic/t/subclass.t
 Basic/t/thread.t
 Basic/t/thread_def.t
+Basic/t/tp-deep.t
+Basic/t/tp-eq_pdl.t
+Basic/t/tp-eq_pdl_clean.t
+Basic/t/tp-import_options.t
+Basic/t/tp-is_pdl.t
 Basic/t/ufunc.t
+Basic/TestDeepPDL.pm
+Basic/TestPDL.pm
 Basic/threads.pm
 Basic/Ufunc/Makefile.PL
 Basic/Ufunc/ufunc.pd
diff --git a/Makefile.PL b/Makefile.PL
index 68bbfab97..fcae79129 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -70,6 +70,7 @@ my %makefile_hash = (
     'IPC::Cmd' => '0.72',
     'Test::Exception' => 0,
     'Test::Warn' => 0,  # for t/pptest.t
+    "Test::Deep" => 0, # for Test::Deep::PDL tests
   },
   BUILD_REQUIRES => {
     'ExtUtils::MakeMaker' => 0,