Herhangi Matematiksel Şeyler: Euler Projesi 261. Soru

Pivot Kare Toplamlar

m > 0 ve n ≥ k tam sayı çifti için k'ye kadar olan (m + 1) ardışık karenin toplamı, (n + 1)'den başlayan m ardışık karenin toplamına eşitse k pozitif tam sayısına bir kare-pivot diyelim: $$(k-m)^2 + ... + k^2 = (n+1)^2 + ... + (n+m)^2.$$
Bazı küçük kare-pivotlar:

4: 3^2 + 4^2 = 5^2
21: 20^2 + 21^2 = 29^2
24: 21^2 + 22^2 + 23^2 + 24^2 = 25^2 + 26^2 + 27^2
110: 108^2 + 109^2 + 110^2 = 133^2 + 134^2

10¹⁰ sayısından küçük tüm farklı kare-pivotların toplamını bulunuz.

Cevap: 238890850232021

Mathematica:
#include<iostream>
#include<set>

using namespace std;

typedef long long Int;

#define REP(i, n) for(int i=0;i<n;i++)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)

const Int MAX=10000000000ll;
//const Int MAX=110;

Int gcd(Int a, Int b){
    if(b==0)return a;
    else return gcd(b, a%b);
}

set<int> st;

int main(){
    Int sum=0;
    for(Int m=1;;m++){
        Int k=m*2*(m+1);
        if(k>MAX)break;
        if(st.find(k)==st.end()){
            st.insert(k);
            sum+=k;
            cout<<k<<" "<<sum<<endl;
        }
    }
    for(Int k1=2;k1*k1*2<=MAX;k1++){
        if(k1%100==0){
            cout<<k1<<endl;
        }
        for(Int n1=k1;;n1++){
            if(gcd(n1, k1)>1)continue;
            Int num=n1*n1-k1*k1, den=k1*k1;
            if(num>den)break;
            for(Int m=1;;m++){
                if(m*num>=den)break;
                Int vd=m*(m+1), vn=(m+1)*k1*k1-m*n1*n1;
                Int h=vd/vn, g=h*(k1+n1), k=g*k1;
                if(k>MAX)break;
                if(vd%vn==0){
                    /*                 Int h=vd/vn;
                                       Int g=h*(k1+n1);
                                       Int k=g*k1;*/
                    if(st.find(k)==st.end()){
                        st.insert(k);
                        sum+=k;
                        cout<<k<<" "<<sum<<endl;
                    }
                }
            }
        }
    }
    cout<<sum<<endl;
    return 0;
}

C/C++:
#include <stdio.h>
#include <math.h>
#include <algorithm>
using namespace std;

typedef __int64 LL;
#define eps 1e-9
#define K 10000000000
#define N 100000
char a[N]={1,1};
int q[N];

LL ar[1000000];
int len=0;

int main()
{
 freopen("out.txt","w",stdout);
 int i,j;
 for(i=1;i<N;i++) q[i]=1;
 for(i=2;i<N;i++) if(a[i]==0)
  for(j=i;j<N;j+=i)
  {
   a[j]=1;
   int k,d=0;
   for(k=j;k%i==0;k/=i) d++;
   for(d=(d+1)/2;d--;) q[j]*=i;
  }
 LL m,y,y_;
 for(m=1;2*m*(m+1)<=K;m++)
 {
  LL M=2*m+1;
  LL maxy=2*m*sqrt(m+1)+eps;
  for(y_=q[m];y_<=maxy;y_+=q[m])
  {
   y=y_;
   LL xx=(m+1)*(y*y/m+1);
   LL x=sqrt(xx)+eps;
   if(x*x==xx && x%2==(m+1)%2 && y%2==m%2)
   {
    for(;;)
    {
     if(x>=y+2*m+1 && (y+m)/2<=K) ar[len++]=(y+m)/2;
     LL z=M*x+(M+1)*y;
     y=(M-1)*x+M*y;
     x=z;
     if((y+m)/2>K) break;
    }
   }
  }
 }
 sort(ar,ar+len);
 LL ans=0;
 for(i=0;i<len;i++)
  if(i==0 || ar[i-1]<ar[i]) ans+=ar[i];
 printf("%I64d\n",ans);
 return 0;
}

ml:
module I = Int64
let (+),( * ),(-),(/) = I.add,I.mul,I.sub,I.div
let etarec2 g = let rec f x y = g f x y in f
let etarec3 g = let rec f x y z = g f x y z in f
let id x = x and (^) f g x = f (g x)
module S = Set.Make(struct type t=int64 let compare = I.compare end)

let truc k l n =
  etarec3 (fun ff n ->  
    let a,b,c,d,e = 2L*n,2L*n+1L,2L*n+2L,2L*n*n-n,2L*n*n in
    let next (x,y) = (x*b+y*a-d,x*c+y*b-e) in
    fun x -> if fst x > l/a then id else
     etarec2 (fun f (x,y) -> if x>l then id else
        (if x>l/a then id else f (next (x,y))) ^ 
        (if y>x+n then ff (y-x) (y-n,y) else id) ^ 
         k x) (next x)) n (n,n)

let answer l =
  let rec f n a = let k,a = truc (fun x (_,y) -> (false,S.add x y))  l n a in 
    if k then (S.fold (+) a 0L) else f (n+1L) (true,a) in
    f 1L (true,S.empty)

let _ = answer 10000000000L

Python:
from math import sqrt
from pr_list import pr  # My list of prime numbers

def f(D):  # returns fundamental solution of X^2 - D Y^2 = 1
    s = int(sqrt(D))
    a1, b1, a2, b2, x, y = 1, 0, 0, 1, 0, 1
    a1, b1, a2, b2, x, y = a1*int((s + x) / y)+a2, b1*int((s + x) / y)+b2, a1, b1, y*int((s + x) / y) - x, (D - (x - y*int((s + x) / y))**2) / y
    while (a1*a1 - D*b1*b1 != 1):
        a1, b1, a2, b2, x, y = a1*int((s + x) / y)+a2, b1*int((s + x) / y)+b2, a1, b1, y*int((s + x) / y) - x, (D - (x - y*int((s + x) / y))**2) / y
    return a1, b1

st = set()

for m in xrange(1, 71000):
    D = m * (m + 1)
    P = 1
    D1 = D
    for i in pr:
        while D1 % (i * i) == 0:
            D1 /= i * i
            P *= i
        if i * i > D1:
            break
    (x0, y0) = f(D1)
    a2, b2 = 2 * m + 1, 2 * D
    while 1:
        assert(b2 * b2 == (a2 * a2 - 1) * m * (m + 1))
        if ((a2 % 2) != 1) or ((b2 % 2) != 0):
            b2, a2 = b2 * x0 + D1 * P * a2 * y0, b2 * y0 + a2 * x0 * P
            assert(a2 % P == 0)
            a2 /= P
            continue
        k = (m * (a2 + 1) + b2) / 2
        if k > 10**10:
            break
        st.add(k)
        b2, a2 = b2 * x0 + D1 * P * a2 * y0, b2 * y0 + a2 * x0 * P
        assert(a2 % P == 0)
        a2 /= P

print sum(st)

Java:
import java.awt.*;
import java.io.*;
import java.math.*;
import java.util.*;
import java.awt.geom.*;
import java.lang.reflect.*;

public class Euler {
 
 public static BigInteger yGivenX(BigInteger m, BigInteger x)
 {
  BigInteger a = m;
  BigInteger b = m.multiply(m.add(BigInteger.ONE));
  BigInteger c = x.multiply(m.multiply(m.add(BigInteger.ONE))).subtract(m.add(BigInteger.ONE).multiply(x).multiply(x));
  
  BigInteger discr = b.multiply(b).subtract(new BigInteger("4").multiply(a).multiply(c));
  
  //need to know if this is a perfect square
  BigInteger lo = BigInteger.ZERO;
  BigInteger hi = discr;
  BigInteger two = new BigInteger("2");
  BigInteger guess = (lo.add(hi)).divide(two);
  while(true)
  {
   int closeness = guess.multiply(guess).compareTo(discr);
   if(closeness == 0)
    break;
   else if(closeness < 0)
   {
    lo = guess;
    BigInteger newGuess = (lo.add(hi)).divide(two);
    if(guess.equals(newGuess))
     return new BigInteger(Long.MAX_VALUE + "");
    else
     guess = newGuess;
   }
   else
   {
    hi = guess;
    BigInteger newGuess = (lo.add(hi)).divide(two);
    if(guess.equals(newGuess))
     return new BigInteger(Long.MAX_VALUE + "");
    else
     guess = newGuess;
   }
  }
  
  return (b.negate().add(guess)).divide(two.multiply(a));
  //(- b +/- guess)/(2a)
  //
 }
 
 public static boolean isSolution(BigInteger m, BigInteger x, BigInteger y)
 {
  BigInteger mPlusOne = m.add(BigInteger.ONE);
  BigInteger mTimesMPlusOne = m.multiply(mPlusOne);
  
  BigInteger first = mPlusOne.multiply(x).multiply(x);
  BigInteger second = m.multiply(y).multiply(y);
  BigInteger third = mTimesMPlusOne.multiply(x);
  BigInteger fourth = mTimesMPlusOne.multiply(y);
  
  return first.subtract(second).subtract(third).subtract(fourth).equals(BigInteger.ZERO);
 }
 
 public static void main(String[] args)
 {  
  BigInteger two = new BigInteger("2");
  
  BigInteger P = new BigInteger("3");
  BigInteger Q = two;
  BigInteger K = BigInteger.ONE;
  BigInteger R = new BigInteger("4");
  BigInteger S = new BigInteger("3");
  
  BigInteger MAX = new BigInteger("10000000000");
  BigInteger sum = BigInteger.ZERO;
  
  TreeSet<Pair<BigInteger,BigInteger>> seeds = new TreeSet<Pair<BigInteger,BigInteger>>();
  HashSet<biginteger> solutions = new HashSet<biginteger>();
  
  java.util.List<long> primes = Primes.listOfPrimesUpTo(100000);
  
  for(BigInteger m = BigInteger.ONE; true; m = m.add(BigInteger.ONE))
  {   
   seeds.clear();
   
   //what are the possible seeds?
   
   //automatic
   seeds.add(new Pair<BigInteger,BigInteger>(BigInteger.ZERO,BigInteger.ZERO));
      
   //x must be greater than y
   
   //if(???)
   {    
    TreeMap<Long,Long> map = new PrimeFactorization(m.intValue(), primes).factorization();
        
    //x must at least contain these primes to become a perfect square
    BigInteger minimal = BigInteger.ONE;
    for(long p: map.keySet())
     minimal = minimal.multiply(new BigInteger(p + ""));
        
    BigInteger xTimes = BigInteger.ONE;
    for(BigInteger x = minimal; x.compareTo(m) < 0; x = minimal.multiply(xTimes))
    {
     BigInteger y = yGivenX(m, x);
     if(y.compareTo(BigInteger.ONE) < 0)
      seeds.add(new Pair<BigInteger,BigInteger>(x,y));
     
     xTimes = xTimes.add(BigInteger.ONE);
    }
    
    
   }
   
   //if(seeds.size() != 1)
   // System.out.println(m + " " + seeds);
   
   for(Pair<BigInteger,BigInteger> seed: seeds)
   {    
    BigInteger xCurrent = seed.getFirst();
    BigInteger yCurrent = seed.getSecond();
    
    while(xCurrent.compareTo(MAX) <= 0)
    {
     if(xCurrent.compareTo(BigInteger.ZERO) > 0  && yCurrent.compareTo(xCurrent) >= 0 && isSolution(m, xCurrent, yCurrent))
     {
      if(solutions.add(xCurrent))
      {
       sum = sum.add(xCurrent);
      }
     }
     
     BigInteger xNew = P.multiply(xCurrent).add(Q.multiply(yCurrent)).add(K);
     BigInteger yNew = R.multiply(xCurrent).add(S.multiply(yCurrent));
     
     xCurrent = xNew;
     yCurrent = yNew;
    }
   }
   
   System.out.println(m + " " + solutions.size() + " " + sum);
   
   
   P = P.add(two);
   Q = Q.add(two);
   K = K.add(BigInteger.ONE);
   R = R.add(two);
   S = S.add(two);
  }

Herhangi Matematiksel Şeyler

Sayfalar

4 Mart 2019

Euler Projesi 261. Soru