Herhangi Matematiksel Şeyler: Euler Projesi 273. Soru

Kareler Toplamı

Şu şekilde tanımlı denklemleri düşünün: $a^2+b^2=N,\, a,b$ ve $N$ tam sayı.

N=65 için iki çözüm mevcut:

$a=1,b=8$ ve $a=4,b=7$.

Yukarıda tanımlanan denklem kümesinin tüm çözümlerinde $a$ değerlerinin toplamına S(N) diyelim.

Bu durumda $S(65)=1+4=5$.

$4k+1<150$ olmak üzere sadece $4k+1$ formundaki asallarla bölünebilen tüm kare-bağımsız $N$ sayıları için $\sum S(N)$ kaçtır?
Cevap: 2032447591196869022

C/C++:
#include<cstdio>
#include<vector>
#include<algorithm>

using namespace std;

typedef long long ll;
typedef pair<ll,ll> pll;
typedef vector<pll> vpll;

const int N = 150;
const int K = 16;
int S[] = {5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149};

pll sum_sqr[N];

ll go(int n, vpll v) {
  if(n == K) {
    ll s = 0;
    for(int i = 0; i < int(v.size()); i++)
      s += min(v[i].first, v[i].second);
    return s;
  }
  vpll new_v;
  int p = S[n];
  for(int i = 0; i < int(v.size()); i++)
    for(int k = -1; k < 2; k += 2)
      {
        pll nv(abs(sum_sqr[p].first * v[i].first +
                   k * sum_sqr[p].second * v[i].second),
               abs(sum_sqr[p].first * v[i].second -
                   k * sum_sqr[p].second * v[i].first));
        if(nv.first > nv.second)
          swap(nv.first, nv.second);
        if(!count(new_v.begin(), new_v.end(), nv))
          new_v.push_back(nv);
      }
  return go(n+1, new_v)+go(n+1, v);
}

int main()
{
  for(int i = 1; i*i < N; i += 2)
    for(int j = 2; i*i+j*j < N; j += 2)
      sum_sqr[i*i+j*j] = pll(i,j);
  printf("%lld\n", go(0, vpll(1, pll(0,1))));
}

ml:
let primes = [
    1.,  2.; 1., 4.; 1., 6.; 1., 10.; 2., 3.; 2., 5.; 2., 7.; 3.,  8.;
    3., 10.; 4., 5.; 4., 9.; 4., 11.; 5., 6.; 5., 8.; 7., 8.; 7., 10.;
]

let main () =
    let rec loop2 acc a b = function
        | []           -> Int64.add acc (Int64.of_float a)
        | (c, d) :: tl ->
            let t0 = b *. d -. a *. c in
            let t1 = b *. c +. a *. d in
            let t2 = b *. c -. a *. d in
            let t3 = a *. c +. b *. d in
            let acc = if t0 < t1 then loop2 acc t0 t1 tl
                                 else loop2 acc t1 t0 tl in
            let acc = loop2 acc (abs_float t2) t3 tl in
                      loop2 acc a b tl in

let rec loop1 acc = function
        | [] -> acc
        | (a, b) :: tl -> loop1 (loop2 acc a b tl) tl in
    loop1 0L primes
;;
Printf.printf "Answer: %Ld\n" (main ())

Matlab:
vecP = primes(150); 
vecP = vecP(mod(vecP,4)==1);
vecA = []; 
vecB = [];
for p = vecP, 
   for a = 1:sqrt(p/2), 
       b = sqrt(p-a^2); 
       if(b == round(b))
          vecA(end+1) = a; 
          vecB(end+1) = b; 
          break;
       end
   end
end
%CELL = cell(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2); 
%CELL(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) = [0;1];
q = zeros(1,length(vecP)); 
sumA = zeros(1,5);
for n = 1:2^length(q)-1,
    q(1) = q(1) + 1; 
    d = 1; 
    while(q(d) == 2) 
       q(d) = 0; 
       q(d+1) = q(d+1)+1; 
       d = d+1;
    end
    disp = polyval(q(end:-1:1),2)
    %p = q+1;
    start = find(q==1,1,'first');
    solnA = vecA(start); 
    solnB = vecB(start);
    for index = (start+1):length(q), 
       if(q(index)) 
          a = vecA(index); 
          b = vecB(index); 
          M = [abs(b*solnA-a*solnB),abs(b*solnA+a*solnB);abs(a*solnA+b*solnB),abs(a*solnA-b*solnB)];
          solnA = min(M);
          solnB = max(M);
       end
    end
    %CELL{p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8),p(9),p(10),p(11),p(12),p(13),p(14),p(15),p(16)} = [solnA;solnB];
    add = zeros(1,5);
    for asoln = solnA, 
       add(1) = add(1) + asoln; 
       for d = 1:4, 
           add(d+1) = add(d+1)+floor(add(d)/10000); 
           add(d) = mod(add(d),10000);
       end
    end
    sumA = sumA + add; 
    for d = 1:4, 
        sumA(d+1) = sumA(d+1)+floor(sumA(d)/10000); 
        sumA(d) = mod(sumA(d),10000);
    end
end

Java:
public class Euler273 {
 final int[] pa = {1,2,1,2,1,4,2,5,3,5,4,1,3,7,4,7};
 final int[] pb = {2,3,4,5,6,5,7,6,8,8,9,10,10,8,11,10};
 // pa^2+pb^2 = primes 5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149

public Euler273(){
  System.out.println( calcS(0,1, 0) );
 }
 private long calcS( long a, long b, int pi ){
  if( pi>=pa.length ) return Math.min(a, b);
  long t = calcS( a, b, pi+1 ); // count all N without factor pi
  long c = pa[pi];
  long d = pb[pi];
  t += calcS( Math.abs(a*d-b*c), a*c+b*d, pi+1);
  if( a!=0 )
   t += calcS( Math.abs(a*c-b*d), a*d+b*c, pi+1);
  return t;
 }
}

Basic:
Sub Euler273()
'About 12 seconds
'It turns out that I stumbled on this by accident - http://en.wikipedia.org/wiki/Brahmagupta-Fibonacci_identity
'p1 = a1^2 + b2^2
'p2 = a2^2 + b2^2
'a and b for p1*p2 become (a1*a2 + b1*b2), abs(a1*b2 - a2*b1)
'                         (a1*b2 + a2*b1), (b1*b2 - a1*a2)
'We want squarefree numbers, so each p can show up only once, giving 65,535 possible numbers

Const MaxNum As Integer = 150
    Dim GoodPrimes() As Double    'first column is the prime, second column is a, third column is b
    Dim a As Long, b As Long, c As Long, Counter
    Dim Bases(1 To 1, 1 To 2) As Double, TotalSum As Variant, StartTime As Single
    
    StartTime = Timer
    ArrayOfPrimesCounter = 0
    Call GetListOfPrimes(MaxNum)
    
    'count how many primes = 1 (mod 4)
    Counter = 0
    For a = 1 To ArrayOfPrimesCounter
        If ArrayOfPrimes(a) Mod 4 = 1 Then
            Counter = Counter + 1
        End If
    Next a
    
    'Find a and b for each prime in the list
    ReDim GoodPrimes(1 To Counter, 1 To 3)
    Counter = 0
    For a = 1 To ArrayOfPrimesCounter
        If ArrayOfPrimes(a) Mod 4 = 1 Then
            Counter = Counter + 1
            GoodPrimes(Counter, 1) = ArrayOfPrimes(a)
            For b = 1 To Sqr(ArrayOfPrimes(a) / 2)
                For c = b + 1 To Sqr(ArrayOfPrimes(a) - b * b)
                    If b * b + c * c = ArrayOfPrimes(a) Then
                        GoodPrimes(Counter, 2) = b
                        GoodPrimes(Counter, 3) = c
                    End If
                Next c
            Next b
        End If
    Next a
    
    TotalSum = 0
    MyCount = 0
    For a = 1 To Counter
        TotalSum = TotalSum + GoodPrimes(a, 2)
        Bases(1, 1) = GoodPrimes(a, 2)
        Bases(1, 2) = GoodPrimes(a, 3)
        
        Call Euler273Recurse(GoodPrimes(a, 1), a + 1, GoodPrimes, Bases, TotalSum)
    Next a
    Debug.Print TotalSum, Timer - StartTime
End Sub
Sub Euler273Recurse(CurNum As Double, CurPos As Integer, GoodPrimes() As Double, CurBases() As Double, ByRef TotalSum As Variant)
    Dim a As Integer, b As Integer, Counter As Long
    Dim a1 As Double, a2 As Double, b1 As Double, b2 As Double, r1 As Double, r2 As Double
    Dim NewBases() As Double

ReDim NewBases(1 To UBound(CurBases, 1) * 2, 1 To 2)        'each new root will double the # of solutions
    For a = CurPos To UBound(GoodPrimes, 1)
        Counter = 0
        For b = 1 To UBound(CurBases, 1)
            Counter = Counter + 1
            a1 = CurBases(b, 1)
            a2 = GoodPrimes(a, 2)
            b1 = CurBases(b, 2)
            b2 = GoodPrimes(a, 3)
            
            'first solution
            r1 = a1 * a2 + b1 * b2
            r2 = Abs(a1 * b2 - a2 * b1)
            If r1 < r2 Then
                NewBases(Counter, 1) = r1
                NewBases(Counter, 2) = r2
            Else
                NewBases(Counter, 1) = r2
                NewBases(Counter, 2) = r1
            End If
            TotalSum = TotalSum + CDec(NewBases(Counter, 1))
            
            'second solution
            Counter = Counter + 1
            r1 = a1 * b2 + a2 * b1
            r2 = b1 * b2 - a1 * a2
            If r1 < r2 Then
                NewBases(Counter, 1) = r1
                NewBases(Counter, 2) = r2
            Else
                NewBases(Counter, 1) = r2
                NewBases(Counter, 2) = r1
            End If
            TotalSum = TotalSum + CDec(NewBases(Counter, 1))
        Next b
        Call Euler273Recurse(CurNum * GoodPrimes(a, 1), a + 1, GoodPrimes, NewBases, TotalSum, MyCount)
    Next a
End Sub

Perl:
#!/usr/bin/perl

use strict;
use warnings;
use Modern::Perl;
use Math::BigInt;
use List::Util qw/ sum /;

my @q = ( 5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137, 149 );
my %ans = (
        5       => [ 1, 2 ],    13      => [ 2, 3 ],    17      => [ 1, 4 ],
        29      => [ 2, 5 ],    37      => [ 1, 6 ],    41      => [ 4, 5 ],
        53      => [ 2, 7 ],    61      => [ 5, 6 ],    73      => [ 3, 8 ],
        89      => [ 5, 8 ],    97      => [ 4, 9 ],    101     => [ 1, 10 ],
        109     => [ 3, 10 ],   113     => [ 7, 8 ],    137     => [ 4, 11 ],
        149     => [ 7, 10 ],
);

my $sum = Math::BigInt->new( 0 );
for ( my $i = 1; $i <= 16; $i++ ) {
        combs( \@q, $i, [ ] );
}
say "Answer to problem 273: $sum";

sub getAnswers {
        my ( $raa, $rb ) = @_;
        return [ $raa ] if scalar @$rb == 0;

my ( $a0, $a1 ) = @$raa;
        my @c = ( );
        for my $ab ( @$rb ) {
                my ( $b0, $b1 ) = @$ab;

my $r = $a0 * $b0 + $a1 * $b1;
                my $s = abs( $a0 * $b1 - $a1 * $b0 );
                push @c, [ $s, $r ];

$r = $a1 * $b1 - $a0 * $b0;
                $s = $a0 * $b1 + $a1 * $b0;
                ( $r, $s ) = ( $s, $r ) if $s < $r;
                push @c, [ $r, $s ];
        }
        return \@c;
}

sub combs {
        my ( $ra, $n, $sol ) = @_;
        my $s;

if ( $n == 0 or scalar @$ra == 0 ) {
                $sum += sum( map { $_->[ 0 ] } @$sol );
                return;
        }
        if ( scalar @$ra == $n ) {
                my $rc = $sol;
                for my $prime ( @$ra ) {
                        $rc = getAnswers( $ans{ $prime }, $rc );
                }

$s = 0;
                for my $pair ( @$rc ) {
                        $s += $pair->[ 0 ];
                }
                $sum += $s;
                return;
        }
        if ( $n == 1 ) {
                $s = 0;
                for my $prime ( @$ra ) {
                        for my $pair ( @{ getAnswers( $ans{ $prime }, $sol ) } ) {
                                $s += $pair->[ 0 ];
                        }
                }
                $sum += $s;
                return;
        }

my @a = @$ra;
        my $i = shift @a;
        combs( \@a, $n - 1, getAnswers( $ans{ $i }, $sol ) );
        combs( \@a, $n, $sol );
}

Pascal:
program project_euler_problem_273;

{ Sum of Squares }

const       z = 150;

var             a : array[1..z] of boolean;
                b : array[1..16,1..2] of int64;
                c : array[1..21523360,1..2] of int64;
            j,k,s : longint;
            g,n,w : int64;

function fn1(x,y,u,v : int64) : int64;
  begin
    fn1 := abs(x*v-y*u);
  end;

function fn2(x,y,u,v : int64) : int64;
  begin
    fn2 := x*u+y*v;
  end;

begin
  a[1] := false;
  for k := 2 to z do
    a[k] := true;
  for k := 2 to trunc(sqrt(z)) do
    if a[k] then
      for s := 2 to z div k do
        a[s*k] := false;
  n := 1;
  for k := 1 to z do
    if a[k] and (k mod 4 = 1) then
      for j := 1 to trunc(sqrt(k/2)) do
        begin
          g := k-j*j;
          s := trunc(sqrt(g));
          if s*s = g then
            begin
              b[n,1] := j;
              b[n,2] := s;
              n := n + 1
            end
        end;
  n := 1;
  c[1,1] := b[1,1];
  c[1,2] := b[1,2];
  for k := 2 to 16 do
    begin
      for j := 1 to n do
        begin
          c[n+j,1] := fn1(c[j,1],c[j,2],b[k,1],b[k,2]);
          c[n+j,2] := fn2(c[j,1],c[j,2],b[k,1],b[k,2]);
          c[2*n+j,1] := fn1(c[j,1],c[j,2],b[k,2],b[k,1]);
          c[2*n+j,2] := fn2(c[j,1],c[j,2],b[k,2],b[k,1])
        end;
      n := 3*n+1;
      c[n,1] := b[k,1];
      c[n,2] := b[k,2]
    end;
  w := 0;
  for k := 1 to n do
    if c[k,1] < c[k,2] then
      w := w + c[k,1]
    else
      w := w + c[k,2];
    write(w);
  readln
end.

Herhangi Matematiksel Şeyler

Sayfalar

25 Eylül 2019

Euler Projesi 273. Soru