prog.pl 1.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768
  1. #!/usr/bin/perl
  2. # Decimal expansion of the sum of reciprocals of Brazilian primes, also called the Brazilian primes constant.
  3. # https://oeis.org/A306759
  4. use 5.014;
  5. #~ use Math::AnyNum qw(:overload float);
  6. #~ my $sum = float(0);
  7. use Math::MPFR;
  8. my $sum = Math::MPFR::Rmpfr_init2(92);
  9. Math::MPFR::Rmpfr_set_ui($sum, 0, 0);
  10. open my $fh, '<', 't5.txt';
  11. my $str = do {
  12. local $/;
  13. <$fh>;
  14. };
  15. #my $data = eval( do {
  16. # local $/;
  17. # <$fh>;
  18. #} );
  19. close $fh;
  20. my $t = Math::MPFR::Rmpfr_init2(92);
  21. while ($str =~ /(\d+)/g) {
  22. #$sum += Math::AnyNum->new($1)->inv;
  23. Math::MPFR::Rmpfr_set_ui($t, $1, 0);
  24. Math::MPFR::Rmpfr_ui_div($t, 1, $t, 0);
  25. Math::MPFR::Rmpfr_add($sum, $sum, $t, 0);
  26. #~ Math::MPFR::Rmpfr_add_d($sum, $sum, 1/$1, 0);
  27. }
  28. say $sum;
  29. # B(10^19) ~ 0.33175446666446018177571079766533
  30. # 10^12 = 0.331754390677227426801521850173028607337819345774
  31. # 10^14 = 0.331754460113696752705452670945987494053288903943
  32. # 10^15 = 0.331754464735448520879667677495078712564872109224
  33. # 10^16 = 0.331754466101486907998310205443082093679034303966
  34. # 10^17 = 0.33175446650734519516960638465
  35. #~ 0.33175446473544852087966761101
  36. #~ 0.33175446473544851633795978262e-1
  37. #~ 0.331754464735448520879667677495078712564872109224
  38. #~ 0.331754464735448516337959901210319518853745198367
  39. #~ 0.33175446473544851633795984906e-1
  40. #~
  41. # Sum of reciprocals of Brazilian primes <~ 10^k, for k=14: 0.33175446011..., k=15: 0.33175446473..., k=16: 0.331754466101..., and k=17: 0.331754466507... (using Maximilian's PARI program from A085104).
  42. __END__
  43. while (<>) {
  44. /^#/ and next;
  45. my ($n, $p) = split(' ');
  46. $p || next;
  47. $p = Math::AnyNum->new($p);
  48. $sum += 1/$p;
  49. say $sum;
  50. }
  51. say $sum;