Euler Projesi 216. Soru

2n²-1 Formundaki Sayıların Asallığı

n>1 olmak üzere t(n)=2n²-1 formundaki t(n) sayılarını ele alalım. Böylesi ilk sayılar 7, 17, 31, 49, 71, 97,  127 ve 161'dir. Burada sadece 49=7*7 ve 161=7*23 asal sayı değildir. n≤10000 için 2202 tane asal olan t(n) sayısı vardır.
n≤50000000 için kaç tane t(n) sayısı vardır?
Cevap: 5437849

Java:

import java.math.*;
public class Problem216 {
 public static void main(String[] args) {
  BigInteger broj=new BigInteger("2");
  BigInteger b2=new BigInteger("2");
  BigInteger b1=new BigInteger("1");
  long koliko=0;
  for (long i=2;i<=50000000;i++){
   if (i%10000==0) System.out.println(i);
   BigInteger b=broj.multiply(broj.multiply(b2));
   if (b.subtract(b1).isProbablePrime(20)) koliko++;
   broj=broj.add(b1);
  }
  System.out.println(koliko);
 }
}

Mathematica:

let np = 50000000
let prime = Array.init (np + 1) (fun i -> 2 * i * i - 1)

let _ =
  let res = ref 0 in
    for i = 2 to np do
      let k = 2 * i * i - 1 in
 if prime.(i) = k then
   incr res;
 let k = prime.(i) in
   if k > 1 then (
     let j = ref (i + k) in
              while !j <= np do
  while prime.(!j) mod k = 0 do
    prime.(!j) <- prime.(!j) / k;
  done;
  j := !j + k
              done;
       let j = ref (-i + k) in
  while !j <= np do
    while prime.(!j) mod k = 0 do
      prime.(!j) <- prime.(!j) / k;
    done;
    j := !j + k
  done
   )
    done;
    Printf.printf "%d\n" !res

 Maple:

/******************************* projectEuler 216   ***************************/
/******************************* Richie 31/10/2012 ****************************/
/******************************* temps :  12.66 s   ***************************/
/******************************* reponse  5437849    **************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include  <math.h>

#define TAILLE_MAX 200000000
#define MAX 50000000
#define qMAX 3*MAX/2
#define llint  long long int

llint projectEuler216(llint max);
llint shanks(llint a, llint q);
llint inverse(llint a,llint q);
llint expomod(llint a,llint n,llint p);
void crible (char *prem, llint taille);

char * prem, *quadratprem;

int main()
{   clock_t temps = clock();
    srand(time(NULL));
    llint i,compteur;

    /******************* creation des tableaux de premiers ****************/

    if ( (prem = malloc(TAILLE_MAX*sizeof(char))) != NULL)
        {  crible(prem,TAILLE_MAX-1);}
        else { exit(1);}
    if ( (quadratprem = malloc((MAX+1)*sizeof(char))) != NULL)
        {
            for (i=0;i<=MAX;i++) quadratprem[i]=1;
        }
        else { exit(1);}

    /******************* fin creation*************************************/

    compteur = projectEuler216(MAX);
    printf("projectEuler216 : %lld\n",compteur);
    printf("temps de calcul : %lf secondes\n", (clock() - temps)/1e6 );

    return 0;
}

/**************** compte les 2n^2 - 1 premiers   **********************/

llint projectEuler216(llint max)
{
    llint n,x0,x1,a,compteur = 0,q;

    for (q=3;q*q<=2*max*max;q++)
        if ( prem[q] )
       {
          a = (q+1)/2;
          if ( expomod(a ,(q-1)/2,q) == 1 )
          {
              x0 = shanks(a,q) ;
              x1 = q - x0;
              if ( q == 2*x0*x0 -1) x0 += q;
              while ( x0 <=MAX)
              {
                  quadratprem[x0] = 0;
                  x0 +=q;
              }
              if ( q == 2*x1*x1 -1) x1 += q;
              while ( x1 <= MAX)
              {
                  quadratprem[x1] = 0;
                  x1 +=q;
              }
          }
       }

       for (n=2;n<=MAX;n++)
          if (quadratprem[n])  compteur++;
        printf("\n");
       return compteur;
}

/**************** calcul d une racine de a modulo q *******************/

llint shanks(llint a, llint q)
{
    llint k=0,r = q-1,u,ksi,rac,alpha,invRac,i,j,a1,b;
    int test = 0;

     while ( r % 2 == 0 )
    {
        r /=2;
        k++;
    }

    while ( !test)
    {
         u = rand () % (q-1) + 1;
         if ( expomod(u, (q-1)/2, q) !=1 ) test = 1;
         ksi = expomod(u,r,q);
    }

    rac =  expomod(a,(r+1)/2,q);
    a = expomod(a,q-1-r,q);
    invRac = 1;
    while (a != 1)
    {
        i = 0;
        a1 = a;
        alpha = 1;
        while ( a1 !=1)
        {
            a1 = ( a1*a1 ) % q;
            alpha =alpha *2;
            i++;
        }
        b = ksi;
        for (j=1; j<k-i; j++)
        {
            b = (b * b ) % q;
        }
        invRac =(invRac * b) % q;
        b = (b * b ) % q;
        a = (a * b ) % q;

    }
    rac = ( rac * inverse(invRac,q) ) %q;
   return rac;
}

/****************** inverse de a modulo b ****************************/

llint inverse(llint a,llint q)
{   llint r0 = a,r1 = q,u0 = 1, v0 = 0, u1 = 0, v1 = 1,r2,u2,v2,quotient;
    while (r1 !=0)
    {
        quotient= r0 / r1;
        r2 = r0 - quotient*r1;
        u2 = u0 - quotient*u1;
        v2 = v0 - quotient*v1;
        r0 = r1; r1 = r2;
        u0 = u1; u1 = u2;
        v0 = v1; v1 = v2;
    }
    while (u0 < 0) u0 +=q;
    return u0;
}

/****************** exponentiation rapide modulo q *******************/

llint expomod(llint a,llint n,llint p)
{
    llint resu = 1;
    while (n != 0)
    {
        if (n % 2 !=0 ) resu = (a *resu ) % p;
        a = (a * a) % p;
        n/=2;
    }
    return resu;
}

/***************************** crible *****************************************/

void crible (char *prem, llint taille)
{
 llint i;
 llint p=2,q,fin;
 for (i=1;i<=taille;i++){ prem[i]=1;}
 prem[0] = 0;prem[1]=0;
 fin = trunc(sqrt(taille));
 while (p<=fin)
 {
  q=p+p;
  while (q<=taille)
  {
   prem[q]=0;
   q=q+p;
  }
  do {p++;} while (!prem[p]);
 }
}

C++:

#include <stdio.h>   
#include <stdlib.h>  
#include <gmp.h>
#include <pthread.h>

void printn(mpz_t n) {
  while(1) {
    gmp_printf("At: %Zd\r", n); fflush(stdout); sleep(1);
  }
}

int main() {
  pthread_t thread;
  pthread_attr_t thread_attr;
  mpz_t n, n2;
  unsigned long result = 0;

  pthread_attr_init(&thread_attr);
  pthread_attr_setdetachstate(&thread_attr, PTHREAD_CREATE_DETACHED);
  
  mpz_init2(n, 1L);
  mpz_init(n2);

  pthread_create(&thread, &thread_attr, (void *(*)(void *)) printn, n);

  while(mpz_cmp_ui(n, 50000000L)<=0) {
    mpz_pow_ui(n2, n, 2UL);
    mpz_mul_ui(n2, n2, 2UL);
    mpz_sub_ui(n2, n2, 1UL);
    if(mpz_probab_prime_p(n2,5)>0) {
      result++;
    }
    mpz_add_ui(n, n, 1);
  }
  printf("\n\nResult: %ul\n", result);
  return 0;
}

Basic:

Private Sub Math_216()

        Dim Number As UInt64
        Dim Base2 As UInt64 = 2
        Dim Base3 As UInt64 = 3
        Dim N As UInt64

        Dim Answer As Int32 = 0

        For N = 2 To 50000000
            Number = 2 * N * N - 1
            If TrialDivide(Number) Or N < 100 Then
                If MR_Primer(Number, Base2) Then
                    If MR_Primer(Number, Base3) Then
                        Answer = Answer + 1
                    End If
                End If
            End If
        Next
        PrintCR(Answer)

    End Sub

 VB.Net:

Public Shared Function Problem216() As String
  Dim count = 0
  Const max = 50 * 1000 * 1000 
  Dim numbers(max) As Long
  For i = 1 To max
    numbers(i) = 2L * i * i - 1
  Next
  For i = 2 To max
    If numbers(i) = 2L * i * i - 1 Then  ' found a prime
      count += 1
    End If
    If numbers(i) > 1 Then  ' found a new factor
      Dim factor = numbers(i)
      For j = i To max Step factor
        Do
          numbers(CInt(j)) \= factor
        Loop While numbers(CInt(j)) Mod factor = 0
      Next
      For j = factor - i To max Step factor
        Do
          numbers(CInt(j)) \= factor
        Loop While numbers(CInt(j)) Mod factor = 0
      Next
    End If
  Next
  Return count.ToString
End Function

Fortran:

 include '\io\!_io.lst'
  include '\timing\!_timing.lst'    !rgl
program p216        ! also see problem 291
    ! This is from st1974's post in the forum for problem 291.
  implicit none
  integer,parameter:: LONG=SELECTED_INT_KIND(16)
  integer(LONG),parameter:: YMAX=5*10_LONG**7
  integer(LONG):: sieve(0:YMAX)
  integer(LONG):: y,p,pp,f
  integer(LONG):: c=0
    real(kind=8) :: time    !rgl
    
    call timeit(time)       !rgl
  sieve(:)=1

  p = 1
  do y=2,YMAX
  !   p=y*y+(y+1)*(y+1)
  !   p = 2*y*(y+1)+1
     p = p + 4*y - 2
     if (sieve(y)==1) then
        c=c+1
        call sieve_factor(y,p)
     else
        pp=p/sieve(y)
        do
            f=gcd(pp,sieve(y))
            if (f==1) exit
            pp=pp/f
        end do
        if (pp>1) call sieve_factor(y,pp)
     end if
  end do

  write(*,*) '>>>',c
  call elapsed(time)        !rgl

 contains 
  
  subroutine sieve_factor(y,p)
    integer(LONG):: y,p,yy
    ! write(*,*) '>',p,y
    ! if (.not. is_prime(p)) stop
    do yy=y,YMAX,p
       sieve(yy)=sieve(yy)*p
    end do
    do yy=p-modulo(y,p),YMAX,p
       sieve(yy)=sieve(yy)*p
    end do
  end subroutine sieve_factor
    

  integer(LONG) function gcd(a,b)
    implicit none
    integer(LONG):: a,b
    integer(LONG):: aa,bb
    
    aa=a
    bb=b
    do
       if (bb==0) then
          gcd=aa
          return
       end if
       aa=modulo(aa,bb)
       if (aa==0) then
          gcd=bb
          return
       end if
       bb=modulo(bb,aa)
    end do
  end function gcd
  
end program p291

!     Output:
!              >>>              5437849
!             elapsed =      9.387 seconds

Delphi:

procedure TForm1.ButtonProblem216Click(Sender: TObject);
const
  Top = 50000000;
var
  n,Count:Integer;
  t:array of Int64;
procedure composite(const p:Int64; u:Integer);
  var v:Int64;
  begin
    v:=u+p;
    while v<=Top do
    begin
      t[v]:=abs(t[v]); while t[v] mod p=0 do t[v]:=t[v] div p;
      inc(v,p);
    end;
    v:=u+p-u shl 1 mod p;
    while v<=Top do
    begin
      t[v]:=abs(t[v]); while t[v] mod p=0 do t[v]:=t[v] div p;
      inc(v,p);
    end;
  end;
begin
  Memo1.Lines.Add(DateTimeToStr(Now()));
  Count:=0;
  setlength(t,Top+1); t[1]:=-1;
  for n:=2 to Top do t[n]:=t[n-1]+2-n shl 2;
  for n:=2 to Top do
  begin
    if t[n]=1 then continue;
    if t[n]<0 then
    begin
      inc(Count); composite(-t[n],n);
    end
    else composite(t[n],n);
  end;
  Memo1.Lines.Add(format('Count(%d) = %d',[Top,Count]));
  Memo1.Lines.Add(DateTimeToStr(Now()));
end;