1#!/usr/bin/env perl
2use strict;
3use warnings;
4
5use Test::More;
6use Math::Prime::Util qw/jordan_totient divisor_sum moebius/;
7
8#my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
9#my $usexs = Math::Prime::Util::prime_get_config->{'xs'};
10#my $usegmp= Math::Prime::Util::prime_get_config->{'gmp'};
11my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
12$use64 = 0 if $use64 && 18446744073709550592 == ~0;
13
14my %jordan_totients = (
15  # A000010
16  1 => [1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12, 10, 22, 8, 20, 12, 18, 12, 28, 8, 30, 16, 20, 16, 24, 12, 36, 18, 24, 16, 40, 12, 42, 20, 24, 22, 46, 16, 42, 20, 32, 24, 52, 18, 40, 24, 36, 28, 58, 16, 60, 30, 36, 32, 48, 20, 66, 32, 44],
17  # A007434
18  2 => [1, 3, 8, 12, 24, 24, 48, 48, 72, 72, 120, 96, 168, 144, 192, 192, 288, 216, 360, 288, 384, 360, 528, 384, 600, 504, 648, 576, 840, 576, 960, 768, 960, 864, 1152, 864, 1368, 1080, 1344, 1152, 1680, 1152, 1848, 1440, 1728, 1584, 2208, 1536],
19  # A059376
20  3 => [1, 7, 26, 56, 124, 182, 342, 448, 702, 868, 1330, 1456, 2196, 2394, 3224, 3584, 4912, 4914, 6858, 6944, 8892, 9310, 12166, 11648, 15500, 15372, 18954, 19152, 24388, 22568, 29790, 28672, 34580, 34384, 42408, 39312, 50652, 48006, 57096],
21  # A059377
22  4 => [1, 15, 80, 240, 624, 1200, 2400, 3840, 6480, 9360, 14640, 19200, 28560, 36000, 49920, 61440, 83520, 97200, 130320, 149760, 192000, 219600, 279840, 307200, 390000, 428400, 524880, 576000, 707280, 748800, 923520, 983040, 1171200],
23  # A059378
24  5 => [1, 31, 242, 992, 3124, 7502, 16806, 31744, 58806, 96844, 161050, 240064, 371292, 520986, 756008, 1015808, 1419856, 1822986, 2476098, 3099008, 4067052, 4992550, 6436342, 7682048, 9762500, 11510052, 14289858, 16671552, 20511148, 23436248, 28629150, 32505856, 38974100, 44015536, 52501944, 58335552, 69343956, 76759038, 89852664, 99168256, 115856200, 126078612, 147008442, 159761600, 183709944, 199526602, 229345006, 245825536, 282458442, 302637500, 343605152, 368321664],
25  # A069091
26  6 => [1, 63, 728, 4032, 15624, 45864, 117648, 258048, 530712, 984312, 1771560, 2935296, 4826808, 7411824, 11374272, 16515072, 24137568, 33434856, 47045880, 62995968, 85647744, 111608280, 148035888, 187858944, 244125000, 304088904, 386889048],
27  # A069092
28  7 => [1, 127, 2186, 16256, 78124, 277622, 823542, 2080768, 4780782, 9921748, 19487170, 35535616, 62748516, 104589834, 170779064, 266338304, 410338672, 607159314, 893871738, 1269983744, 1800262812, 2474870590, 3404825446],
29);
30
31my @A001615 = (1,3,4,6,6,12,8,12,12,18,12,24,14,24,24,24,18,36,20,36,32,36,24,48,30,42,36,48,30,72,32,48,48,54,48,72,38,60,56,72,42,96,44,72,72,72,48,96,56,90,72,84,54,108,72,96,80,90,60,144,62,96,96,96,84,144,68,108,96);
32
33plan tests => scalar(keys %jordan_totients)
34              + 2  # Dedekind psi calculated two ways
35              + 2  # Calculate J5 two different ways
36              + 2 * $use64 # Jordan totient example
37              ;
38
39###### Jordan Totient
40while (my($k, $tref) = each (%jordan_totients)) {
41  my @tlist = map { jordan_totient(0+$k, $_) } 1 .. scalar @$tref;
42  is_deeply( \@tlist, $tref, "Jordan's Totient J_$k" );
43}
44
45{
46  my @psi_viaj;
47  my @psi_viamobius;
48  foreach my $n (1 .. scalar @A001615) {
49    push @psi_viaj, int(jordan_totient(2, $n) / jordan_totient(1, $n));
50    push @psi_viamobius, int($n * divisor_sum( $n, sub { moebius($_[0])**2 / $_[0] } ) + 0.5);
51  }
52  is_deeply( \@psi_viaj, \@A001615, "Dedekind psi(n) = J_2(n)/J_1(n)" );
53  is_deeply( \@psi_viamobius, \@A001615, "Dedekind psi(n) = divisor_sum(n, moebius(d)^2 / d)" );
54}
55
56{
57  my $J5 = $jordan_totients{5};
58  my @J5_jordan = map { jordan_totient(5, $_) } 1 .. scalar @$J5;
59  is_deeply( \@J5_jordan, $J5, "Jordan totient 5, using jordan_totient");
60  my @J5_moebius = map { my $n = $_; divisor_sum($n, sub { my $d=shift; $d**5 * moebius($n/$d); }) } 1 .. scalar @$J5;
61  is_deeply( \@J5_moebius, $J5, "Jordan totient 5, using divisor sum" );
62}
63
64if ($use64) {
65  is( jordan_totient(4, 12345), 22902026746060800, "J_4(12345)" );
66  # Apostal page 48, 17a.
67  is( divisor_sum( 12345, sub { jordan_totient(4,$_[0]) } ),
68      # was int(12345 ** 4), but Perl 5.8.2 gets it wrong.
69      int(12345*12345*12345*12345),
70      "n=12345, k=4  :   n**k = divisor_sum(n, jordan_totient(k, d))" );
71}
72