Herhangi Matematiksel Şeyler: Euler Projesi 258. Soru

Gecikmeli Fibonacci Dizisi

Bu dizi aşağıdaki şekilde tanımlanıyor:

g_k = 1, 0 ≤ k ≤ 1999 ise
g_k = g_k-2000 + g_k-1999, k ≥ 2000 ise.

k = 10¹⁸ için g_k mod 20092010 değerini bulunuz.
Cevap: 12747994

C/C++:
#include <cstdio>

#include <vector>

using namespace std;

typedef long long int llint;

// const static llint K = 300000;
const static llint K = 1000000000000000000LL;
const static int M = 20092010;
const static int N = 2000;

typedef vector<llint> TIntVector;
typedef vector<tintvector> TIntMatrix;

int Stupido(llint k) {
    const static int SIZE = 2*N;
    TIntVector v;
    v.resize(SIZE);
    for (size_t i = 0; i < N; ++i)
        v[i] = 1;
    for (size_t i = N; i <= k; ++i)
        v[i % SIZE] = (v[(i - N + 1) % SIZE] + v[(i - N) % SIZE]) % M;
    return v[k % SIZE];
}

void MakeMatrix(int n, int m, TIntMatrix* result) {
    TIntVector dummy(m, 0);
    result->resize(n, dummy);
}

void MulMatrix(const TIntMatrix& a, const TIntMatrix& b, TIntMatrix* c) {
    MakeMatrix(a.size(), b.size(), c);
    for (size_t i = 0; i < a.size(); ++i) {
        for (size_t j = 0; j < a[0].size(); ++j) {
            if (a[i][j]) {
                for (size_t k = 0; k < b.size(); ++k) {
                    (*c)[i][k] += a[i][j]*b[j][k];
                    (*c)[i][k] %= M;
                }
            }
        }
    }
}

void MakeMatrixId(int m, TIntMatrix* result) {
    MakeMatrix(m, m, result);
    for (size_t i = 0; i < m; ++i)
        (*result)[i][i] = 1;
}

void PowerMatrix(const TIntMatrix& m, llint power, TIntMatrix* result) {
    TIntMatrix res;
    MakeMatrixId(m.size(), &res);
    for (llint i = 0; i < power; ++i) {
        TIntMatrix tmp;
        MulMatrix(res, m, &tmp);
        res = tmp;
    }
    *result = res;
}

void FastPowerMatrix(const TIntMatrix& m, llint power, TIntMatrix* result) {
    TIntMatrix a;
    MakeMatrixId(m.size(), &a);
    TIntMatrix b = m;
    while (power) {
        printf("power: %lld\n", power);
        if (power & 1) {
            TIntMatrix tmp;
            MulMatrix(a, b, &tmp);
            a = tmp;
            --power;
        } else {
            TIntMatrix tmp;
            MulMatrix(b, b, &tmp);
            b = tmp;
            power >>= 1;
        }
    }
    *result = a;
}

int Fast(llint k) {
    TIntMatrix r;
    MakeMatrix(N, N, &r);
    for (int i = 0; i < N - 1; ++i)
        r[i][i + 1] = 1;
    r[N - 1][0] = 1;
    r[N - 1][1] = 1;
    
    FastPowerMatrix(r, k - N + 1, &r);
    
    TIntMatrix v;
    MakeMatrix(N, 1, &v);
    for (int i = 0; i < N; ++i)
        v[i][0] = 1;
    TIntMatrix res;
    MulMatrix(r, v, &res);
    
    return res[N - 1][0];
}

int main() {
    // printf("%d\n", Stupido(K));
    printf("%d\n", Fast(K));

return 0;
}

Python:
#include <math.h>
#include <stdio.h>

long long int A1[2000*2000], A[2000*2000], C[2000*2000];

void printmat(long long int *M, int P) {

int i,j;
  for (i=0; i<P; i++) {
    for (j=0; j<P; j++) {
      printf("%4lld ",M[i*P+j]);
    }
    printf("\n");
  }
  return;
}

int main() {

long long int tmp,n,N;
  int k,kmax,P,M,i,j,m;
  int B[100];

n=1000000000000000000LL;
  P=2000; 
  M=20092010;
  //  P=5; 
  //  M=13;

printf("P=%d, M=%d, n=%lld\n",P,M,n);

for (i=0;i<2000;i++) {
    for (j=0;j<2000;j++) {
      A[P*i+j]=0LL;
      A1[P*i+j]=0LL;
      C[P*i+j]=0LL;
    }
  }

for (i=0; i<P; i++) C[i*P+i]=1LL;
  
  for (i=0; i<P; i++) A1[i*P+i]=1LL;
  for (i=1; i<P; i++) A1[i*P+i-1]=1LL;
  A1[P-1]=1LL;
  A1[2*P-1]=1LL;

//  N=(n+P)/P;
  N=n/P;
  kmax=100;
  for (k=63; k>=0; k--) {
    if (pow(2,k)<=N) {
      if (kmax==100) kmax=k;
      B[k]=1;
      N-=pow(2,k);
    } else {
      B[k]=0;
    }
  }

//  N=n+P;
  N=n;
  printf("kmax: %d  N: %lld \n",kmax,N);
  for (k=0; k<=kmax; k++) printf("%d ",B[k]);
  printf("\n");
  
  for (k=0; k<=kmax; k++) {
    printf("k=%d B[%d]=%d\n",k,k,B[k]);
    if (B[k]==1) {
      //A=C*A1
      for (i=0; i<P; i++) {
        for (j=0; j<P; j++) {
          A[i*P+j]=0;
          for (m=0; m<P; m++) {
     A[i*P+j]=(A[i*P+j]+(C[i*P+m]*A1[m*P+j]))%M;
          }
 }
      }
      //C=A
      for (i=0; i<P; i++)
        for (j=0; j<P; j++)
          C[i*P+j]=A[i*P+j];
    }
    // A=A1*A1
    for (i=0; i<P; i++) {
      for (j=0; j<P; j++) {
 A[i*P+j]=0;
 for (m=0; m<P; m++) {
   A[i*P+j]=(A[i*P+j]+(A1[i*P+m]*A1[m*P+j]))%M;
 }
      }
    }
    // A1=A
    for (i=0; i<P; i++)
      for (j=0; j<P; j++)
 A1[i*P+j]=A[i*P+j];
    
  }
  //  printf("C=\n");
  //  printmat(C,P);
  tmp=0LL;
  for (i=0; i<P; i++)
    tmp=(tmp+C[i*P+N%P])%M;
  
  printf(" Answer: %lld %lld \n",N,tmp);
}

Haskell:
p :: Integer
p = 20092010

add :: Integer -> Integer -> Integer
add x y = (x + y) `mod` p

shift :: (Num a) => Int -> [a] -> [a]
shift k xs = let (hd, tl) = splitAt k $ reverse xs
             in zipWith (+) (reverse hd ++ reverse tl) (0:reverse hd ++ repeat 0)

mul :: [Integer] -> [Integer] -> [Integer]
mul as bs = map (`mod` p).foldl1 add' $ [shift k $ map (b*) as | (b, k) <- zip bs [0..]]
    where add' xs ys = strict $ zipWith add xs ys
          strict xs = foldr seq xs xs

pow :: (Integral t) => [Integer] -> t -> [Integer]
pow as 1 = as
pow as n | even n    = map (`mod` p) $ pow (mul as as) (div n 2)
         | otherwise = mul as $ pow as (n-1)

main :: IO ()
main = print.foldl1 add'.pow (1:1:replicate 1998 0) $ 5*10^14
    where add' x y | x + y >= p = x+y-p
                   | otherwise  = x+y

Java:
import static java.lang.Math.*;
import static java.lang.System.*;

/**
 * Euler Problem 258
 * @author kevinsogo
 */
public class P258 {
    static long pow(long b, int e) {//fast exponentiation (for 10^18).
        if (e == 0) return 1;
        if ((e & 1) == 0) return pow(b * b, e >> 1);
        return b * pow(b, e - 1);
    }
    public static void main (String args[]) throws Exception {
        long time = nanoTime();
        long k = args.length == 0 ? pow(10, 18) : Long.parseLong(args[0]);
        final int n = 2000;
        final long mod = 20092010;
        
        long b[] = new long[n];//b contains the n base cases.
        long c[] = new long[n << 1];//c is a temp array to contain new base cases.
        long p[] = new long[n + 1];//p represents the coefficients of the current polynomial. Initially 1 - x^1999 - x^2000
        long q[] = new long[n + 1];//q is a temp array to contain the new p.
        
        // initialize stuff
        java.util.Arrays.fill(b, 1);
        p[0] = 1;
        p[1999] = p[2000] = -1;
        
        // do this while k is not in one of our base cases.
        for (int m = n << 1; k >= n; k >>= 1) {
        
            // just to keep track of time
            out.printf("get(%d): %.8f s.\n", k, (nanoTime() - time) * 1e-9); 
            
            // copy the first n values.
            for (int i = 0; i < n; i++) {
                c[i] = b[i];
            }
            
            // compute the next n values using the recurrence (derived from coeffs of p).
            for (int i = n; i < m; i++) {
                c[i] = 0;
                for (int j = 1; j <= n; j++) {
                    if (p[j] == 0) continue;//skip (p[j] = 0 has no effect).
                    c[i] -= p[j] * c[i - j];
                }
                c[i] %= mod;
            }
            // c contains the new base values.
            //copy the even/odd terms as the new base cases (depending on the parity of k).
            for (int i = 0, j = (int)(k & 1); i < n; i++, j += 2) {
                b[i] = c[j];
            }
            
            // set q to be the product of p and it's "conjugate".
            // we store only even-power terms since odd-power terms are all zero.
            for (int i = 0; i <= n; i++) {
                q[i] = p[i] * p[i];
                int lim = min(n - i, i);
                for (int j = 1, a = -2; j <= lim; j++, a = -a) {
                    q[i] += p[i + j] * p[i - j] * a;
                }
                if ((i & 1) != (n & 1)) {
                    q[i] *= -1;
                }
            }
            
            //take q as the new polynomial.
            for (int i = 1; i <= n; i++) {
                p[i] = q[i] % mod;
            }
        }
        
        //now that k < n, the answer is now in the base cases.
        long get = b[(int)k] % mod;
        get = (get + mod) % mod;
        out.printf("%d in %.8f s.\n", get, (nanoTime() - time) * 1e-9);
    }
}

ml:
module I = Int64
let ( *:),(+:),(%:),(/:) = I.mul, I.add, I.rem, I.div

let table = Array.init 3999 (fun i ->
  if i<2000 then [i] else if i<3999 then [(i-2000); (i-1999)] else [0;1;1999])

let ((+:),( *:)) = let m = 20092010L in (fun a b -> (a+:b) %: m),(fun a b -> (a*:b)%:m)

let ( ** ) a b =
 let c = Array.make 2000 0L in
 for i=0 to 1999 do if a.(i)>0L then
  for j=0 to 1999 do if b.(j)>0L then
   let d = a.(i)*:b.(j) in List.iter (fun k -> c.(k)<-c.(k)+:d) table.(i+j) done done; c
   
let exp = let rec f a b c = if c=0L then a else f (if c%:2L = 1L then a**b else a) (b**b) (c/:2L) in
 f (Array.init 2000 (fun i -> if i=0 then 1L else 0L))

let nth n = Array.fold_left (+:) 0L (exp (Array.init 2000 (fun i -> if i=1 then 1L else 0L)) n)

let answer = nth 1000000000000000000L

Fortran:
include "Restklassen.f90"

program p258
  use Restklassen ! module containing the Chinese Remainder
                  ! function RSynchron
  implicit none
  integer,parameter:: LONG=SELECTED_INT_KIND(16)
  integer,parameter:: NDIM=2000
  integer:: p

type Mvector  
     integer:: c(0:NDIM-1)  ! coefficients of linear combination
  end type Mvector
  type(Mvector):: M_pow_np(0:NDIM-1)  ! precompute M^(n*p) mod p

integer(LONG):: k=10_LONG**18
  type(Restklasse):: r2,r5,r859,r2339,r

call set_p(2)
  r2=RK(g(k),2)
  write(*,*) r2

call set_p(5)
  r5=RK(g(k),5)
  write(*,*) r5

call set_p(859)
  r859=RK(g(k),859)
  write(*,*) r859

call set_p(2339)
  r2339=RK(g(k),2339)
  write(*,*) r2339

r=RSynchron(r2,r5)
  r=RSynchron(r,r859)
  r=RSynchron(r,r2339)
  write(*,*) r

contains

integer function g(k) ! g(k) mod p
    implicit none
    integer(LONG):: k,kk
    type(Mvector):: MM,Mprod
    integer:: n,i,r,r1,i1

Mprod=M_pow_n(0)
    kk=k
    n=0
    r1=0
    do while (kk>0)
       r=mod(kk,p+0_LONG)

if (r>0) then
          if (r/=r1) then    ! otherwise reuse MM
             MM=M_pow_n(r)
             i1=0
          end if
          do i=i1+1,n
             MM=pow_p(MM)
          end do
          Mprod=mul(MM,Mprod)
          r1=r
          i1=i
       end if

kk=kk/p
       n=n+1
    end do
    g = mod( sum(Mprod%c(:)) , p )

end function g

subroutine inc(z,x)
    implicit none
    integer:: z,x
    if (x==0) return
    z=mod(z+x,p)
  end subroutine inc

type(Mvector) function mul(x,y) result(z) ! z=x*y
    implicit none
    type(Mvector):: x,y
    integer:: i,j
    z%c(:)=0
    do j=0,NDIM-1
       if (y%c(j)==0) cycle ! speed up if y is sparse
       do i=0,NDIM-1-j
          call inc(z%c(i+j) ,x%c(i)*y%c(j))
       end do
    end do
    do j=1,NDIM-1
       if (y%c(j)==0) cycle ! speed up if y is sparse
       do i=NDIM-j,NDIM-1
          call inc(z%c(i+j-NDIM)  , x%c(i)*y%c(j))
          call inc(z%c(i+j-NDIM+1), x%c(i)*y%c(j))
       end do
    end do
  end function mul

type(Mvector) function M_pow_n(n) result(z) ! z=M^n (for small n)
    implicit none
    integer:: n
    z%c(:)=0
    if (n<NDIM) then
       z%c(n)=1
    else if (n<2*NDIM-1) then
       z%c(n-NDIM)=1
       z%c(n-NDIM+1)=1
    else
       write(*,*) 'M_pow_n out of range'
       stop
    end if
  end function M_pow_n

type(Mvector) function pow_p(x) result(z) ! z=x^p
    implicit none
    type(Mvector):: x
    integer:: i,n
    z%c(:)=0
    do n=0,NDIM-1
       if (x%c(n)==0) cycle
       do i=0,NDIM-1
          call inc(z%c(i), x%c(n) * M_pow_np(n)%c(i))
       end do
    end do
  end function pow_p

subroutine set_p(p0)
    implicit none
    integer:: p0,n
    p=p0
    M_pow_np(0) = M_pow_n(0)
    M_pow_np(1) = M_pow_n(p)
    do n=2,NDIM-1
       M_pow_np(n) = mul(M_pow_np(n-1) , M_pow_np(1))
    end do 
  end subroutine set_p

end program p258

Herhangi Matematiksel Şeyler

Sayfalar

2 Şubat 2019

Euler Projesi 258. Soru