From aa2b425a9dae2dac26eaa298d3b782a1c822123b Mon Sep 17 00:00:00 2001 From: Alex Jordan Date: Mon, 25 May 2026 22:01:56 -0700 Subject: [PATCH 1/2] change the Round function to combat machine rounding error --- macros/core/PGauxiliaryFunctions.pl | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/macros/core/PGauxiliaryFunctions.pl b/macros/core/PGauxiliaryFunctions.pl index b46cb0443f..c06dbed667 100644 --- a/macros/core/PGauxiliaryFunctions.pl +++ b/macros/core/PGauxiliaryFunctions.pl @@ -175,12 +175,29 @@ =head2 Round Round(1.789,2) returns 1.79 +This macro assumes an input x having lots of consecutive 9s such as 1.4999999999999997 is already +the result of machine rounding error, and was meant to be 1.5. And so will round up to 2, not down +to 1. Do not use this macro if a problem might work with honest numbers like 1.4999999999999997 +and need to round it. + =cut -# Round contributed bt Mark Schmitt 3-6-03 sub Round { - if (@_ == 1) { $_[0] > 0 ? int $_[0] + 0.5 : int $_[0] - 0.5 } - elsif (@_ == 2) { $_[0] > 0 ? Round($_[0] * 10**$_[1]) / 10**$_[1] : Round($_[0] * 10**$_[1]) / 10**$_[1] } + my ($x, $n) = @_; + $n = 0 unless $n; + my $e = (split(/E/, sprintf("%E", $x)))[1] + 0; # exponent for $x + my $s = ($x < 0 ? -1 : 1); # the sign of $x + my $N = $e + $n; # number of digits to retain + $x += $s * 10**($e - 15); # adjust for repeated 9s + return sprintf("%.${N}E", $x) + 0 unless $N < 0; # round the adjusted value + + # For zero original digits, we add a digit just above the first one + # in $x, round that, then remove the added digit, getting 0 if $x + # didn't round up, or 1 in the proper place if it did. This means + # that 0.005 rounds to .01, for example, when 2 digits are requested. + + my $d = $s * 10**($e + 1); + return sprintf("%.0E", $x + $d) - $d; } =head2 lcm From a29ec6898325e6463615c05ba5dbccb3c5462041 Mon Sep 17 00:00:00 2001 From: Peter Staab Date: Mon, 1 Jun 2026 09:44:28 -0400 Subject: [PATCH 2/2] added test for rounding and added subtests. --- t/macros/pgaux.t | 207 ++++++++++++++++++++++++----------------------- 1 file changed, 108 insertions(+), 99 deletions(-) diff --git a/t/macros/pgaux.t b/t/macros/pgaux.t index f6977ea347..3573f28963 100644 --- a/t/macros/pgaux.t +++ b/t/macros/pgaux.t @@ -13,105 +13,114 @@ do "$ENV{PG_ROOT}/t/build_PG_envir.pl"; loadMacros('PGauxiliaryFunctions.pl'); -# test step functions - -is(step(8), 1, "step: positive number"); -is(step(-8), 0, "step: negative number"); -is(step(0), 0, "step: step(0)=0"); - -# test floor function - -is(floor(0.5), 0, "floor: positive non-integer"); -is(floor(-0.5), -1, "floor: negative non-integer"); -is(floor(1), 1, "floor: positive integer"); -is(floor(0), 0, "floor: floor(0)=0"); -is(floor(-1), -1, "floor: negative integer"); - -# test ceiling function - -is(ceil(0.5), 1, "ceil: positive non-integer"); -is(ceil(-0.5), 0, "ceil: negative non-integer"); -is(ceil(1), 1, "ceil: positive integer"); -is(ceil(0), 0, "ceil: floor(0)=0"); -is(ceil(-1), -1, "ceil: negative integer"); - -# max/min functions - -is(max(1, 2, 3, 9, 4, 5, 6, 8), 9, "max: set of integers"); -is(max(0.1, -2.3, 1.345, 2.71712, -1000.1), 2.71712, "max: set of decimals"); -is(min(1, 2, 3, 9, 4, 5, 6, 8), 1, "min: set of integers"); -is(min(0.1, -2.3, 1.345, 2.71712, -1000.1), -1000.1, "min: set of decimals"); - -# round function - -is(round(0.95), 1, "round: fractional part > 0.5"); -is(round(0.45), 0, "round: fractional part < 0.5"); -is(round(0.5), 1, "round: fractional part = 0.5"); -is(round(-0.95), -1, "round: fractional part > 0.5 and negative"); -is(round(-0.45), 0, "round: fractional part < 0.5 and negative"); -is(round(-0.5), -1, "round: fractional part = 0.5 and negative"); - -# Round function which takes a second number, the number of digits to round to - -is(Round(1.793, 2), 1.79, "Round to 2 digits: test 1"); -is(Round(1.797, 2), 1.80, "Round to 2 digits: test 2"); -is(Round(1.795, 2), 1.80, "Round to 2 digits: test 3"); -is(Round(-1.793, 2), -1.79, "Round to 2 digits: test 1"); -is(Round(-1.797, 2), -1.80, "Round to 2 digits: test 2"); -is(Round(-1.795, 2), -1.80, "Round to 2 digits: test 3"); - -is(Round(15.793, -1), 20, "Round to -1 digits (nearest 10)"); - -# lcm - -is(lcm(20, 30), 60, "lcm: non relatively prime numbers"); -is(lcm(5, 6), 30, "lcm: relatively prime numbers"); -is(lcm(2, 3, 4), 12, "lcm: 3 numbers"); -is(lcm(2, 3, 4, 5, 6, 7, 8), 840, "lcm: 7 numbers"); - -# gcd -is(gcd(16, 8), 8, "gcd: 2 powers of 2"); -is(gcd(10, 9), 1, "gcd: 2 relatively prime"); - -is(gcd(10, 20, 30, 40), 10, "gcd: 4 multiples of 10"); - -# isPrime -is(isPrime(7), 1, "isPrime: 7 is prime"); -is(isPrime(2), 1, "isPrime: 2 is prime"); -is(isPrime(15), 0, "isPrime: 15 is not prime"); +subtest 'Step function' => sub { + is(step(8), 1, "step: positive number"); + is(step(-8), 0, "step: negative number"); + is(step(0), 0, "step: step(0)=0"); +}; + +subtest 'Floor function' => sub { + is(floor(0.5), 0, "floor: positive non-integer"); + is(floor(-0.5), -1, "floor: negative non-integer"); + is(floor(1), 1, "floor: positive integer"); + is(floor(0), 0, "floor: floor(0)=0"); + is(floor(-1), -1, "floor: negative integer"); +}; + +subtest 'Ceiling function' => sub { + is(ceil(0.5), 1, "ceil: positive non-integer"); + is(ceil(-0.5), 0, "ceil: negative non-integer"); + is(ceil(1), 1, "ceil: positive integer"); + is(ceil(0), 0, "ceil: floor(0)=0"); + is(ceil(-1), -1, "ceil: negative integer"); +}; + +subtest 'Min and Max functions' => sub { + is(max(1, 2, 3, 9, 4, 5, 6, 8), 9, "max: set of integers"); + is(max(0.1, -2.3, 1.345, 2.71712, -1000.1), 2.71712, "max: set of decimals"); + is(min(1, 2, 3, 9, 4, 5, 6, 8), 1, "min: set of integers"); + is(min(0.1, -2.3, 1.345, 2.71712, -1000.1), -1000.1, "min: set of decimals"); +}; + +subtest 'round and Round functions' => sub { + is(round(0.95), 1, "round: fractional part > 0.5"); + is(round(0.45), 0, "round: fractional part < 0.5"); + is(round(0.5), 1, "round: fractional part = 0.5"); + is(round(-0.95), -1, "round: fractional part > 0.5 and negative"); + is(round(-0.45), 0, "round: fractional part < 0.5 and negative"); + is(round(-0.5), -1, "round: fractional part = 0.5 and negative"); + + # Round function which takes a second number, the number of digits to round to + + is(Round(1.793, 2), 1.79, "Round to 2 digits: test 1"); + is(Round(1.797, 2), 1.80, "Round to 2 digits: test 2"); + is(Round(1.795, 2), 1.80, "Round to 2 digits: test 3"); + is(Round(-1.793, 2), -1.79, "Round to 2 digits: test 1"); + is(Round(-1.797, 2), -1.80, "Round to 2 digits: test 2"); + is(Round(-1.795, 2), -1.80, "Round to 2 digits: test 3"); + + is(Round(15.793, -1), 20, "Round to -1 digits (nearest 10)"); + + # tests that round handles some cases related to precision loss in calculations + + is(Round(134.49999999999997, 0), 135.0, 'Round a number close to 0.5'); + is(Round(0.01499999999991, 3), 0.015, 'Round a number close to 0.005 to 2 digits'); +}; + +subtest 'lcm and gcd functions' => sub { + is(lcm(20, 30), 60, "lcm: non relatively prime numbers"); + is(lcm(5, 6), 30, "lcm: relatively prime numbers"); + is(lcm(2, 3, 4), 12, "lcm: 3 numbers"); + is(lcm(2, 3, 4, 5, 6, 7, 8), 840, "lcm: 7 numbers"); + + # gcd + is(gcd(16, 8), 8, "gcd: 2 powers of 2"); + is(gcd(10, 9), 1, "gcd: 2 relatively prime"); + + is(gcd(10, 20, 30, 40), 10, "gcd: 4 multiples of 10"); + +}; + +subtest 'isPrime function' => sub { + is(isPrime(7), 1, "isPrime: 7 is prime"); + is(isPrime(2), 1, "isPrime: 2 is prime"); + is(isPrime(15), 0, "isPrime: 15 is not prime"); +}; # random_coprime - -my $sum = 0; -for my $i (1 .. 100) { - my @coprimes = random_coprime([ 1 .. 20 ], [ 1 .. 20 ]); - $sum += gcd($coprimes[0], $coprimes[1]); -} -is($sum, 100, "random_coprime: 100 tests in 1..20,1..20"); - -$sum = 0; - -for my $i (1 .. 100) { - my @coprimes = random_coprime([ -9 .. -1, 1 .. 9 ], [ 1 .. 9 ], [ 1 .. 9 ]); - $sum += gcd(@coprimes); -} -is($sum, 100, "random_coprime: 100 tests in [-9..-1,1..9],[1..9],[1..9]"); - -my ($sum1, $sum2, $sum3, $sum4) = (0, 0, 0); -for my $i (1 .. 100) { - my @coprimes = random_pairwise_coprime([ -9 .. -1, 1 .. 9 ], [ 1 .. 9 ], [ 1 .. 9 ]); - $sum1 += gcd(@coprimes); - $sum2 += gcd($coprimes[0], $coprimes[1]); - $sum3 += gcd($coprimes[0], $coprimes[2]); - $sum4 += gcd($coprimes[1], $coprimes[2]); -} -is($sum1 + $sum2 + $sum3 + $sum4, 400, "random_pairwise_coprime: 100 tests of [-9..-1,1..9],[1..9],[1..9]"); - -# reduce -# it would be nicer to directly compare the arrays -my @my_arr = (3, 4); -my @res = reduce(15, 20); -is($my_arr[0], $res[0], "reduce: correct numerator"); -is($my_arr[1], $res[1], "reduce: correct denominator"); - +subtest 'random_coprime function' => sub { + my $sum = 0; + for my $i (1 .. 100) { + my @coprimes = random_coprime([ 1 .. 20 ], [ 1 .. 20 ]); + $sum += gcd($coprimes[0], $coprimes[1]); + } + is($sum, 100, "random_coprime: 100 tests in 1..20,1..20"); + + $sum = 0; + + for my $i (1 .. 100) { + my @coprimes = random_coprime([ -9 .. -1, 1 .. 9 ], [ 1 .. 9 ], [ 1 .. 9 ]); + $sum += gcd(@coprimes); + } + is($sum, 100, "random_coprime: 100 tests in [-9..-1,1..9],[1..9],[1..9]"); + + my ($sum1, $sum2, $sum3, $sum4) = (0, 0, 0); + for my $i (1 .. 100) { + my @coprimes = random_pairwise_coprime([ -9 .. -1, 1 .. 9 ], [ 1 .. 9 ], [ 1 .. 9 ]); + $sum1 += gcd(@coprimes); + $sum2 += gcd($coprimes[0], $coprimes[1]); + $sum3 += gcd($coprimes[0], $coprimes[2]); + $sum4 += gcd($coprimes[1], $coprimes[2]); + } + is($sum1 + $sum2 + $sum3 + $sum4, 400, "random_pairwise_coprime: 100 tests of [-9..-1,1..9],[1..9],[1..9]"); +}; + +subtest 'reduce function' => sub { + + # Note: this is testing reduction of fractions. Not sure why this is here. + my @my_arr = (3, 4); + my @res = reduce(15, 20); + is($my_arr[0], $res[0], "reduce: correct numerator"); + is($my_arr[1], $res[1], "reduce: correct denominator"); +}; done_testing;