Lucas-Lehmer Primality Test

ll_test.cpp

// Lucas-Lehmer Primality Test

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <stdint.h>

const int modulus=30; // example modulus=6,12,18,30

bool ll_test(uint64_t P)
{
    // ll test for P>2*modulus return true if 2^P-1 is prime
    uint64_t bW=(uint64_t)modulus;
    int nR;
    std::vector<int> RW; //vectors of the remainders associated with possible classes using (modulo bW)
    int i_P=0;
    if (P>2*bW)
    {
        int pt=2;
        int pmax;
        int d=1;        
        std::vector<char> Remainder_t(modulus,true); 
        while (pt*pt<=modulus)
        {
            if (modulus%pt==0 && Remainder_t[pt]==true)
            {
                pmax=pt;
                for (int j=pt*pt;j< modulus;j+=pt)
                    Remainder_t[j]=false;
            }
            pt+=d;
            d=2;
        }
        RW.push_back(1);
        for (int j=pmax+1;j< modulus;j++)
            if (Remainder_t[j]==true)
                RW.push_back(j);
        nR=RW.size();

        while(P%bW!=RW[i_P] && i_P<nR)
            i_P++;
    }
    else
        std::cout<<std::endl<<"Enter P> "<< 2*bW<<std::endl;
	
    if (P>2*bW && i_P<nR)
    {    
        uint64_t bW_mod=(1<<modulus);
        std::vector<uint64_t> mod_f(nR,0);
        std::vector<uint64_t> coef_f(nR,0);
        for (int i=0;i<nR;i++)
        {
            mod_f[i]=1<<RW[i]; 
            coef_f[i]=1<<RW[nR-1-i]; 
        }            

        uint64_t dim_v=(P+(uint64_t)RW[nR-1-i_P])/bW;
        std::vector<uint64_t> prod_t(dim_v,0);
        std::vector<uint64_t> esp_L_i_m2(dim_v+1,0);
        std::vector<uint64_t> coef_L_i_m2(dim_v,0);


        uint64_t l_m2=1;      
        esp_L_i_m2[0]=0;
        coef_L_i_m2[0]=2;     //L_1=S_0-2=2

        uint64_t    coef_sum,coef_t,esp_t,c_t,i1,j1;

        for (uint64_t step=2;step<P;step++)
        {
            //begin multiplication and partial operation (modulo 2^P-1) - L_i+1 =L_i *L_i +4 *L_i
            for (i1=0;i1<l_m2 && 2*esp_L_i_m2[i1]<dim_v;i1++)
            {
                prod_t[esp_L_i_m2[i1]]+=4*coef_L_i_m2[i1];
                c_t=1;
                for (j1=i1;j1<l_m2 && (esp_L_i_m2[i1]+esp_L_i_m2[j1])<dim_v-1;j1++)
                {
                    esp_t=esp_L_i_m2[i1]+esp_L_i_m2[j1];
                    coef_t=(c_t*coef_L_i_m2[i1]*coef_L_i_m2[j1]);
                    prod_t[esp_t]+=coef_t%bW_mod;
                    prod_t[esp_t+1]+=coef_t/bW_mod;
                    c_t=2;
                }
                if (j1<l_m2 &&  (esp_L_i_m2[i1]+esp_L_i_m2[j1])==dim_v-1)
                {
                    coef_t=(c_t*coef_L_i_m2[i1]*coef_L_i_m2[j1]);
                    prod_t[0]+=coef_t/mod_f[i_P];
                    prod_t[dim_v-1]+= coef_t%mod_f[i_P];
                    c_t=2;
                    j1++;
                }
                for (;j1<l_m2;j1++)
                {
                    esp_t=esp_L_i_m2[i1]+esp_L_i_m2[j1]-dim_v;
                    coef_t=(c_t*coef_L_i_m2[i1]*coef_L_i_m2[j1]);
                    prod_t[esp_t]+=(coef_f[i_P]*(coef_t%bW_mod))%bW_mod;
                    prod_t[esp_t+1]+=coef_f[i_P]*(coef_t/bW_mod)+(coef_f[i_P]*(coef_t%bW_mod))/bW_mod;
                    c_t=2;
                }
            }
            for (;i1<l_m2;i1++)
            {
                prod_t[esp_L_i_m2[i1]]+=4*coef_L_i_m2[i1];
                c_t=1;
                for (j1=i1;j1<l_m2;j1++)
                {
                    esp_t=esp_L_i_m2[i1]+esp_L_i_m2[j1]-dim_v;
                    coef_t=(c_t*coef_L_i_m2[i1]*coef_L_i_m2[j1]);
                    prod_t[esp_t]+=(coef_f[i_P]*(coef_t%bW_mod))%bW_mod;
                    prod_t[esp_t+1]+=coef_f[i_P]*(coef_t/bW_mod)+(coef_f[i_P]*(coef_t%bW_mod))/bW_mod;
                    c_t=2;
                }
            }
            //end multiplication and partial operation (modulo 2^P-1) 

            //begin completion of the (modulo 2^P-1) operation 
            do
            {
                l_m2=0;
                coef_sum=0;
                prod_t[0]+=prod_t[dim_v-1]/mod_f[i_P];
                prod_t[dim_v-1]%= mod_f[i_P];
                for (uint64_t i=1;i<dim_v;i++)
                {
                    prod_t[i]+=prod_t[i-1]/bW_mod;
                    prod_t[i-1]%=bW_mod;
                    if (prod_t[i-1]>0)
                    {
                        esp_L_i_m2[l_m2]=(i-1);
                        coef_L_i_m2[l_m2]=prod_t[i-1];
                        coef_sum+=prod_t[i-1];
                        l_m2++;
                    }
                }
                if (prod_t[dim_v-1]>0 && prod_t[dim_v-1]<mod_f[i_P])
                {
                    esp_L_i_m2[l_m2]=dim_v-1;
                    coef_L_i_m2[l_m2]=prod_t[dim_v-1];
                    l_m2++;
                }
            }
            while (prod_t[dim_v-1]>=mod_f[i_P]);
           //end completion of the (modulo 2^P-1) operation
 
            std::fill(prod_t.begin(), prod_t.end(), 0);
        }
        if (l_m2==dim_v)
        {
            if ( coef_sum==((bW_mod-1)*dim_v-bW_mod-1) && coef_L_i_m2[0]==bW_mod-3 && coef_L_i_m2[dim_v-1]==mod_f[i_P]-1)
               return true;
        }
    }
    return false;
}

int main()
{ 
    // ll_test(P) -- ll test for P>2*modulus 
    std::cout <<"ll test (1 prime, 0 composite): "<< ll_test(607)<< std::endl; 

    return 0;
}