123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150 |
- #!/usr/bin/perl
- # a(n) is the least prime p for which the continued fraction expansion of sqrt(p) has exactly n consecutive 1's starting at position 2.
- # https://oeis.org/A307453
- # a(30) = 326757090818617
- # Algorithm from:
- # http://web.math.princeton.edu/mathlab/jr02fall/Periodicity/mariusjp.pdf
- use 5.010;
- use strict;
- use warnings;
- use ntheory qw(is_prime);
- use Math::GMPz;
- sub isok {
- my ($n, $want) = @_;
- $n = Math::GMPz->new("$n");
- my $x = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_sqrt($x, $n);
- return 0 if Math::GMPz::Rmpz_perfect_square_p($n);
- my $y = Math::GMPz::Rmpz_init_set($x);
- my $z = Math::GMPz::Rmpz_init_set_ui(1);
- my $r = Math::GMPz::Rmpz_init();
- my @cfrac = ($x);
- Math::GMPz::Rmpz_add($r, $x, $x); # r = x+x
- my $count = 0;
- do {
- my $t = Math::GMPz::Rmpz_init();
- # y = (r*z - y)
- Math::GMPz::Rmpz_submul($y, $r, $z); # y = y - t*z
- Math::GMPz::Rmpz_neg($y, $y); # y = -y
- # z = floor((n - y*y) / z)
- Math::GMPz::Rmpz_mul($t, $y, $y); # t = y*y
- Math::GMPz::Rmpz_sub($t, $n, $t); # t = n-t
- Math::GMPz::Rmpz_divexact($z, $t, $z); # z = t/z
- # t = floor((x + y) / z)
- Math::GMPz::Rmpz_add($t, $x, $y); # t = x+y
- Math::GMPz::Rmpz_tdiv_q($t, $t, $z); # t = floor(t/z)
- $r = $t;
- #push @cfrac, $t;
- ++$count;
- if ($count <= $want) {
- Math::GMPz::Rmpz_cmp_ui($t, 1) == 0 or return 0;
- }
- elsif ($count == $want+1) {
- return (Math::GMPz::Rmpz_cmp_ui($t, 1) != 0);
- }
- } until (Math::GMPz::Rmpz_cmp_ui($z, 1) == 0);
- return 0;
- }
- #~ use ntheory qw(forprimes);
- #~ my $k = 1;
- #~ my $want = 2;
- #~ local $| = 1;
- #~ while (1) {
- #~ ++$k;
- #~ if (isok($k, $want)) {
- #~ print "$k, ";
- #~ ++$want;
- #~ $k = 1;
- #~ }
- #~ }
- #~ __END__
- # Limit_{n->infinity} (sqrt(a(n)) - floor(sqrt(a(n)))) = A094214.
- use Math::AnyNum qw(phi round floor PREC 1024);
- my $k = 1;
- my $want = 3;
- local $| = 1;
- while (1) {
- ++$k;
- my $n = round(($k + phi - 1)**2);
- # is_prime($n) || next;
- if (isok($n, $want)) {
- print "$n, ";
- $k = 1;
- ++$want;
- }
- }
- # 7, 13, 244, 58, 135, 3921, 819, 2081, 68444, 13834, 35955, 1220181, 244647, 639389, 21861404, 4374866, 11448871, 392143681, 78439683
- # 7, 13, 3797, 5273, 4987, 90371, 79873, 2081, 111301, 1258027, 5325101, 12564317, 9477889, 47370431, 709669249, 1529640443, 2196104969, 392143681, 8216809361, 30739072339, 200758317433,
- # 7, 13, 3797, 5273, 4987, 90371, 79873, 2081, 111301, 1258027, 5325101, 12564317, 9477889, 47370431, 709669249, 1529640443, 2196104969, 392143681
- # 7, 13, 3797, 5273, 4987, 90371, 79873, 2081, 111301, 1258027, 5325101, 12564317, 9477889, 47370431, 709669249, 1529640443, 2196104969, 392143681
- #my @arr = ( 3, 31, 7, 13, 3797, 5273, 4987, 90371, 79873, 2081, 111301, 1258027, 5325101, 12564317, 9477889, 47370431, 709669249, 1529640443, 2196104969, 392143681);
- #foreach my $k(0..$#arr) {
- # say isok($arr[$k], $k+1);
- #}
- __END__
- 7, 13, 3797, 5273, 4987, 90371, 79873, 2081, 111301, 1258027, 5325101, 12564317, 9477889, 47370431, 709669249, 1529640443, 2196104969, 392143681, 8216809361, 30739072339, 200758317433, 370949963971, 161356959383, 1788677860531, 7049166342469, 4484287435283, 3690992602753, ^C
- perl r.pl 286.67s user 0.16s system 99% cpu 4:48.23 total
- sqrt(1) = [1]
- sqrt(2) = [1, 2]
- sqrt(3) = [1, 1, 2]
- sqrt(4) = [2]
- sqrt(5) = [2, 4]
- sqrt(6) = [2, 2, 4]
- sqrt(7) = [2, 1, 1, 1, 4]
- sqrt(8) = [2, 1, 4]
- sqrt(9) = [3]
- sqrt(10) = [3, 6]
- sqrt(11) = [3, 3, 6]
- sqrt(12) = [3, 2, 6]
- sqrt(13) = [3, 1, 1, 1, 1, 6]
- sqrt(14) = [3, 1, 2, 1, 6]
- sqrt(15) = [3, 1, 6]
- sqrt(16) = [4]
- sqrt(17) = [4, 8]
- sqrt(18) = [4, 4, 8]
- sqrt(19) = [4, 2, 1, 3, 1, 2, 8]
- sqrt(20) = [4, 2, 8]
|