#!/usr/bin/perl -w use strict; use Math::BigInt try => 'GMP'; my \$DBL_DECIMAL_DIG= 17; my \$resultprecision= \$DBL_DECIMAL_DIG; if( !@ARGV ) { print < computes mod 2 pi accurately EOF exit; } # from http://mapage.noos.fr/echolalie/l127.htm : my \$pistr= <new(sprintf("%f", \$num)); # turn bit shifts into powers of 10 \$gmpnum->blsft(\$numshift, 5); my \$dbl2pi= 6.28318530717958647692; my \$quot= \$num / \$dbl2pi; # The numerical error in the modulo operation below is \$quot * ULP(\$twopi), # so we add an appropriate number of decimal digits of additional accuracy. my \$accdecimals= int(log(\$quot) / log(10)) + 1; # >= 1 \$gmpnum->blsft(\$resultprecision + \$accdecimals, 10); my \$pidecimals= 1 + \$resultprecision + \$accdecimals + \$numshift; \$pistr= substr(\$pistr, 0, \$pidecimals); my \$twopi= Math::BigInt->new("\$pistr"); \$twopi->bmul(2); \$gmpnum->bmod(\$twopi); # add .5*ULP to round up and strip surplus digits \$gmpnum->brsft(\$accdecimals + \$numshift - 1, 10); \$gmpnum->badd(5); \$gmpnum->brsft(1, 10); my \$resultstr= \$gmpnum->bstr(); if( length(\$resultstr) > \$resultprecision ) { print substr(\$resultstr, 0, length(\$resultstr) - \$resultprecision), ".", substr(\$resultstr, -\$resultprecision), "\n"; } else { print "0."; \$resultstr= "0" x (\$resultprecision - length(\$resultstr)) . \$resultstr; print substr(\$resultstr, 0, \$resultprecision), "\n"; }