#!/usr/bin/perl -w

use strict;

use Math::BigInt try => 'GMP';

my $DBL_DECIMAL_DIG= 17;
my $resultprecision= $DBL_DECIMAL_DIG;

if( !@ARGV ) {
    print <<EOF;
usage: mod2pi.pl <number>
computes <number> mod 2 pi accurately
EOF
    exit;
}

# from http://mapage.noos.fr/echolalie/l127.htm :
my $pistr= <<EOF;
3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510
58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 82148
08651 32823 06647 09384 46095 50582 23172 53594 08128 48111 74502
84102 70193 85211 05559 64462 29489 54930 38196 44288 10975 66593
34461 28475 64823 37867 83165 27120 19091 45648 56692 34603 48610
45432 66482 13393 60726 02491 41273 72458 70066 06315 58817 48815
20920 96282 92540 91715 36436 78925 90360 01133 05305 48820 46652
13841 46951 94151 16094 33057 27036 57595 91953 09218 61173 81932
61179 31051 18548 07446 23799 62749 56735 18857 52724 89122 79381
83011 94912 98336 73362 44065 66430 86021 39494 63952 24737 19070
21798 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132
00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 14684
40901 22495 34301 46549 58537 10507 92279 68925 89235 42019 95611
21290 21960 86403 44181 59813 62977 47713 09960 51870 72113 49999
99837 29780 49951 05973 17328 16096 31859 50244 59455 34690 83026
42522 30825 33446 85035 26193 11881 71010 00313 78387 52886 58753
32083 81420 61717 76691 47303 59825 34904 28755 46873 11595 62863
88235 37875 93751 95778 18577 80532 17122 68066 13001 92787 66111
95909 21642 01989 38095 25720 10654 85863 27886 59361 53381 82796
82303 01952 03530 18529 68995 77362 25994 13891 24972 17752 83479
13151 55748 57242 45415 06959 50829 53311 68617 27855 88907 50983
81754 63746 49393 19255 06040 09277 01671 13900 98488 24012 85836
16035 63707 66010 47101 81942 95559 61989 46767 83744 94482 55379
77472 68471 04047 53464 62080 46684 25906 94912 93313 67702 89891
52104 75216 20569 66024 05803 81501 93511 25338 24300 35587 64024
74964 73263 91419 92726 04269 92279 67823 54781 63600 93417 21641
21992 45863 15030 28618 29745 55706 74983 85054 94588 58692 69956
90927 21079 75093 02955 32116 53449 87202 75596 02364 80665 49911
98818 34797 75356 63698 07426 54252 78625 51818 41757 46728 90977
77279 38000 81647 06001 61452 49192 17321 72147 72350 14144 19735
68548 16136 11573 52552 13347 57418 49468 43852 33239 07394 14333
45477 62416 86251 89835 69485 56209 92192 22184 27255 02542 56887
67179 04946 01653 46680 49886 27232 79178 60857 84383 82796 79766
81454 10095 38837 86360 95068 00642 25125 20511 73929 84896 08412
84886 26945 60424 19652 85022 21066 11863 06744 27862 20391 94945
04712 37137 86960 95636 43719 17287 46776 46575 73962 41389 08658
32645 99581 33904 78027 59009 94657 64078 95126 94683 98352 59570
98258 22620 52248 94077 26719 47826 84826 01476 99090 26401 36394
43745 53050 68203 49625 24517 49399 65143 14298 09190 65925 09372
21696 46151 57098 58387 41059 78859 59772 97549 89301 61753 92846
81382 68683 86894 27741 55991 85592 52459 53959 43104 99725 24680
84598 72736 44695 84865 38367 36222 62609 91246 08051 24388 43904
51244 13654 97627 80797 71569 14359 97700 12961 60894 41694 86855
58484 06353 42207 22258 28488 64815 84560 28506 01684 27394 52267
46767 88952 52138 52254 99546 66727 82398 64565 96116 35488 62305
77456 49803 55936 34568 17432 41125 15076 06947 94510 96596 09402
52288 79710 89314 56691 36867 22874 89405 60101 50330 86179 28680
92087 47609 17824 93858 90097 14909 67598 52613 65549 78189 31297
84821 68299 89487 22658 80485 75640 14270 47755 51323 79641 45152
37462 34364 54285 84447 95265 86782 10511 41354 73573 95231 13427
16610 21359 69536 23144 29524 84937 18711 01457 65403 59027 99344
03742 00731 05785 39062 19838 74478 08478 48968 33214 45713 86875
19435 06430 21845 31910 48481 00537 06146 80674 91927 81911 97939
95206 14196 63428 75444 06437 45123 71819 21799 98391 01591 95618
14675 14269 12397 48940 90718 64942 31961 56794 52080 95146 55022
52316 03881 93014 20937 62137 85595 66389 37787 08303 90697 92077
34672 21825 62599 66150 14215 03068 03844 77345 49202 60541 46659
25201 49744 28507 32518 66600 21324 34088 19071 04863 31734 64965
14539 05796 26856 10055 08106 65879 69981 63574 73638 40525 71459
10371 20628 04390 39759 51567 71577 00420 33786 99360 07230 55876
31763 59421 87312 51471 20532 92819 18261 86125 86732 15791 98414
84882 91644 70609 57527 06957 22091 75671 16722 91098 16909 15280
17350 67127 48583 22287 18352 09353 96572 51210 83579 15136 98820
91444 21006 75103 34671 10314 12671 11369 90865 85163 98315 01970
16515 11685 17143 76576 18351 55650 88490 99898 59982 38734 55283
31635 50764 79185 35893 22618 54896 32132 93308 98570 64204 67525
90709 15481 41654 98594 61637 18027 09819 94309 92448 89575 71282
89059 23233 26097 29971 20844 33573 26548 93823 91193 25974 63667
30583 60414 28138 83032 03824 90375 89852 43744 17029 13276 56180
93773 44403 07074 69211 20191 30203 30380 19762 11011 00449 29321
51608 42444 85963 76698 38952 28684 78312 35526 58213 14495 76857
26243 34418 93039 68642 62434 10773 22697 80280 73189 15441 10104
46823 25271 62010 52652 27211 16603 96665 57309 25471 10557 85376
34668 20653 10989 65269 18620 56476 93125 70586 35662 01855 81007
29360 65987 64861 17910 45334 88503 46113 65768 67532 49441 66803
96265 79787 71855 60845 52965 41266 54085 30614 34443 18586 76975
14566 14068 00700 23787 76591 34401 71274 94704 20562 23053 89945
61314 07112 70004 07854 73326 99390 81454 66464 58807 97270 82668
30634 32858 78569 83052 35808 93306 57574 06795 45716 37752 54202
11495 57615 81400 25012 62285 94130 21647 15509 79259 23099 07965
47376 12551 76567 51357 51782 96664 54779 17450 11299 61489 03046
39947 13296 21073 40437 51895 73596 14589 01938 97131 11790 42978
28564 75032 03198 69151 40287 08085 99048 01094 12147 22131 79476
47772 62241 42548 54540 33215 71853 06142 28813 75850 43063 32175
18297 98662 23717 21591 60771 66925 47487 38986 65494 94501 14654
06284 33663 93790 03976 92656 72146 38530 67360 96571 20918 07638
32716 64162 74888 80078 69256 02902 28472 10403 17211 86082 04190
00422 96617 11963 77921 33757 51149 59501 56604 96318 62947 26547
36425 23081 77036 75159 06735 02350 72835 40567 04038 67435 13622
22477 15891 50495 30984 44893 33096 34087 80769 32599 39780 54193
41447 37744 18426 31298 60809 98886 87413 26047 21569 51623 96586
45730 21631 59819 31951 67353 81297 41677 29478 67242 29246 54366
80098 06769 28238 28068 99640 04824 35403 70141 63149 65897 94092
43237 89690 70697 79422 36250 82216 88957 38379 86230 01593 77647
16512 28935 78601 58816 17557 82973 52334 46042 81512 62720 37343
14653 19777 74160 31990 66554 18763 97929 33441 95215 41341 89948
54447 34567 38316 24993 41913 18148 09277 77103 86387 73431 77207
54565 45322 07770 92120 19051 66096 28049 09263 60197 59882 81613
32316 66365 28619 32668 63360 62735 67630 35447 76280 35045 07772
35547 10585 95487 02790 81435 62401 45171 80624 64362 67945 61275
31813 40783 30336 25423 27839 44975 38243 72058 35311 47711 99260
63813 34677 68796 95970 30983 39130 77109 87040 85913 37464 14428
22772 63465 94704 74587 84778 72019 27715 28073 17679 07707 15721
34447 30605 70073 34924 36931 13835 04931 63128 40425 12192 56517
98069 41135 28013 14701 30478 16437 88518 52909 28545 20116 58393
41965 62134 91434 15956 25865 86557 05526 90496 52098 58033 85072
24264 82939 72858 47831 63057 77756 06888 76446 24824 68579 26039
53527 73480 30480 29005 87607 58251 04747 09164 39613 62676 04492
56274 20420 83208 56611 90625 45433 72131 53595 84506 87724 60290
16187 66795 24061 63425 22577 19542 91629 91930 64553 77991 40373
40432 87526 28889 63995 87947 57291 74642 63574 55254 07909 14513
57111 36941 09119 39325 19107 60208 25202 61879 85318 87705 84297
25916 77813 14969 90090 19211 69717 37278 47684 72686 08490 03377
02424 29165 13005 00516 83233 64350 38951 70298 93922 33451 72201
38128 06965 01178 44087 45196 01212 28599 37162 31301 71144 48464
09038 90644 95444 00619 86907 54851 60263
EOF
$pistr =~ tr/0-9//cd;

my $num= 0.0 + $ARGV[0];
if( $num < 0 ) {
    print "- ";
    $num= -$num;
}
if( $num <= 6.28 ) {
    print "$num\n";
    exit;
}

sprintf("%e", $num) =~ /e\+(\d+)$/ or die "Could not match exponential notation";
my $numdecimals= $1 + 1;
my $numshift= 0;
if( $numdecimals < $DBL_DECIMAL_DIG ) {
    # We want to convert to a large integer type in order to perform an
    # integer modulo operation, but numbers as small as this may have
    # fractional digits.  So we multiply them (exactly) with a power of 2 to
    # obtain an integer.
    $numshift= int((($DBL_DECIMAL_DIG - $numdecimals) * 7 + 1) / 2);
    # 2 ** 3.5 = 11.3...
    $num *= 2 ** $numshift;
    # exact multiplication resulting in integer
}
# bug in constructor means we need the sprintf()
my $gmpnum= Math::BigInt->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";
}

