123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 24 October 2017
- # https://github.com/trizen
- # Algorithm for finding solutions to the equation a^2 + b^2 = n,
- # for any given positive integer `n` for which such a solution exists.
- # The number of returned solutions is at least as many as
- # the number of unique prime factors p = 1 (mod 4) in `n`.
- # For numbers with primes powers p^k = 1 (mod 4), for k > 1, not all the possible solutions are returned.
- # For example, when n = 9925 = 5^2 * 397, only the following two solutions are returned: [58, 81], [33, 94].
- # The missing solution for 9925, is: [30, 95].
- # This algorithm is efficient when the factorization of `n` is known.
- # See also:
- # https://oeis.org/A001481
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use ntheory qw(sqrtmod factor_exp chinese forsetproduct);
- sub sum_of_two_squares_solution ($n) {
- $n == 0 and return [0, 0];
- my $prod1 = 1;
- my $prod2 = 1;
- my @prime_powers;
- foreach my $f (factor_exp($n)) {
- if ($f->[0] % 4 == 3) { # p = 3 (mod 4)
- $f->[1] % 2 == 0 or return; # power must be even
- $prod2 *= $f->[0]**($f->[1] >> 1);
- }
- elsif ($f->[0] == 2) { # p = 2
- if ($f->[1] % 2 == 0) { # power is even
- $prod2 *= $f->[0]**($f->[1] >> 1);
- }
- else { # power is odd
- $prod1 *= $f->[0];
- $prod2 *= $f->[0]**(($f->[1] - 1) >> 1);
- push @prime_powers, $f->[0];
- }
- }
- else { # p = 1 (mod 4)
- $prod1 *= $f->[0]**$f->[1];
- push @prime_powers, $f->[0]**$f->[1];
- }
- }
- $prod1 == 1 and return [$prod2, 0];
- $prod1 == 2 and return [$prod2, $prod2];
- my %table;
- foreach my $pp (@prime_powers) {
- my $r = sqrtmod($pp - 1, $pp);
- push @{$table{$pp}}, [$r, $pp], [$pp - $r, $pp];
- }
- my @square_roots;
- forsetproduct {
- push @square_roots, chinese(@_);
- } values %table;
- my @solutions;
- foreach my $r (@square_roots) {
- my $s = $r;
- my $q = $prod1;
- while ($s * $s > $prod1) {
- ($s, $q) = ($q % $s, $s);
- }
- push @solutions, [$prod2 * $s, $prod2 * ($q % $s)];
- }
- return sort { $a->[0] <=> $b->[0] } do {
- my %seen;
- grep { !$seen{$_->[0]}++ } map {
- [sort { $a <=> $b } @$_]
- } @solutions;
- };
- }
- foreach my $n (1 .. 1e5) {
- (my @solutions = sum_of_two_squares_solution($n)) || next;
- say "$n = " . join(' = ', map { "$_->[0]^2 + $_->[1]^2" } @solutions);
- # Verify solutions
- foreach my $solution (@solutions) {
- if ($n != $solution->[0]**2 + $solution->[1]**2) {
- die "error for $n: (@$solution)\n";
- }
- }
- }
- __END__
- 999826 = 99^2 + 995^2 = 315^2 + 949^2 = 699^2 + 715^2 = 525^2 + 851^2
- 999828 = 318^2 + 948^2
- 999844 = 410^2 + 912^2 = 312^2 + 950^2
- 999848 = 62^2 + 998^2
- 999850 = 43^2 + 999^2 = 321^2 + 947^2
- 999853 = 387^2 + 922^2
- 999857 = 544^2 + 839^2 = 401^2 + 916^2
- 999860 = 698^2 + 716^2 = 154^2 + 988^2
- 999869 = 262^2 + 965^2 = 613^2 + 790^2
- 999881 = 484^2 + 875^2 = 341^2 + 940^2
- 999882 = 309^2 + 951^2 = 651^2 + 759^2
- 999890 = 421^2 + 907^2 = 473^2 + 881^2
- 999892 = 324^2 + 946^2
- 999898 = 697^2 + 717^2 = 213^2 + 977^2
- 999909 = 678^2 + 735^2 = 222^2 + 975^2
- 999914 = 667^2 + 745^2
- 999917 = 109^2 + 994^2
- 999937 = 89^2 + 996^2 = 44^2 + 999^2
- 999938 = 77^2 + 997^2
- 999940 = 696^2 + 718^2 = 126^2 + 992^2 = 448^2 + 894^2 = 178^2 + 984^2
- 999941 = 446^2 + 895^2 = 370^2 + 929^2
- 999944 = 638^2 + 770^2
- 999946 = 585^2 + 811^2
- 999949 = 243^2 + 970^2 = 450^2 + 893^2
- 999952 = 444^2 + 896^2
- 999953 = 568^2 + 823^2
- 999954 = 375^2 + 927^2 = 327^2 + 945^2
- 999956 = 500^2 + 866^2
- 999961 = 644^2 + 765^2
- 999962 = 541^2 + 841^2 = 239^2 + 971^2
- 999968 = 452^2 + 892^2
- 999970 = 627^2 + 779^2 = 247^2 + 969^2
- 999973 = 658^2 + 753^2 = 118^2 + 993^2 = 63^2 + 998^2 = 622^2 + 783^2
- 999981 = 141^2 + 990^2
- 999986 = 365^2 + 931^2 = 695^2 + 719^2
- 999997 = 194^2 + 981^2 = 454^2 + 891^2
- 1000000 = 352^2 + 936^2
|