#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <limits.h>

int primes[] =
{
    2, 3, 5, 7, 11, 13, 17, 19,
    23, 29, 31, 37, 41, 43, 47, 53,
    59, 61, 67, 71, 73, 79, 83, 89,
    97, 101, 103, 107, 109, 113, 127, 131,
    137, 139, 149, 151, 157, 163, 167, 173,
    179, 181, 191, 193, 197, 199, 211, 223,
    227, 229, 233, 239, 241, 251, 257, 263,
    269, 271, 277, 281, 283, 293, 307, 311,INT_MAX
};
#define PI 3.1415926535897932384626433832795

// 计算x以内y光滑数的个数
// 参考文献
// 1. 《On integers free of large prime factors》 by A.Hildebrand and G.Tenenbaum. .
// 2. 《A Fast Algorithm for Approximately Counting Smooth Numbers》 by Jonathan P. Sorenson

/*
  用到的几个函数
  1. zeta(s,y)
   // prod= (1- p1^(-s))(1- p2^(-s))...(1- pi^(-s)), pi<=y
   // return 1/prod
  2. phi1(s,y)= sum(log(p)/( p^s-1)), p<=y
  3. phi2(s,y)= sum(p^s*(log p)^2/( p^s-1)^2), p<=y
  4. 函数 alpha(s,c,v),表示方程v - phi1(s,c) =0 的解s
  5. 函数 HT(x,y,s)= (x^s * zeta(s,y))/(s * sqrt(pi * phi2(s,y))
  */

  //1. zeta(s, y)
  // prod= (1- p1^(-s))(1- p2^(-s))...(1- pi^(-s)), pi<=y
  // return 1/prod
double zeta(double s, int y)
{
    double prod = 1.0;
    for (int i = 0; primes[i] <= y; i++)
    {
        double p = (double)primes[i];
        double item = 1.0 - pow(p, -s);
        prod *= item;
    }
    return 1.0 / prod;
}

// 2. phi1(s, y) = sum( log(p)/(p^s - 1) ) for p<=y
double phi1(double s, int y)
{
    double sum = 0.0;
    for (int i = 0; primes[i] <= y; i++)
    {
        double p = (double)primes[i];
        double item = log(p) / (pow(p, s) - 1.0);
        sum += item;
    }
    return sum;
}

//3. phi2(s, y) = sum( p^s * (log p)^2 / (p^s - 1)^2 ) for p <= y
double phi2(double s, int y)
{
    double sum = 0.0;
    for (int i = 0; primes[i] <= y; i++)
    {
        double p = (double)primes[i];
 double ps1= pow(p,s)-1;
        double item = pow(p,s)*log(p)*log(p)/(ps1*ps1);
        sum += item;
    }
    return sum;
}

typedef struct alpha_iter_st
{
    double a0;
    double v0;
    double a1;
    double v1;
    double v; // target value
}ALPHA_ST;


/*
 函数 phi1(x,c),x为自变量，c是常数，
 v0=phi1(x0,c),v1=phi1(x1,c),
 求 方程 v - phi1(x, c) = 0 的解x
 假定函数 phi1(x,c)在区间[x0,x1]，近似为一次函数， 则
  斜率 k=(v1-v0)/(x1-x0)
=> (x-x1)=k/(v-v1)
=> x= x1+(x-x1)
=> x= x1+(v-v0)/k
*/

// 计算函数 phi1(s,y)的近似值s0,
// s0= log(1 + y/v)/log(y), 这里y是素数
inline double get_s0_v1(double v, int y)
{
    return log(1.0 + (double)y/v)/log((double)y);
}

/*
   首先用牛顿迭代法求 函数4.
   方法。
   1.计算 初值 r0
   2.使用牛顿迭代法求r1,
   3.若r1-r0<1e-12, 结束计算，返回r1
      否则 r0=r1, 继续重复步骤2
*/

//inline double alpha_core(ALPHA_ST* p)
// return max(p),p是素数，p<=y,
int correct(int y)
{
    if (y < 2)
        return -1;
    
    int i = 0;
    while (primes[i] <= y)
        i++;
    return primes[i - 1];
}

double get_x1(int y, double v, double x0, double v0, double* pv)
{
    double v1;
    double x1 = x0;

    if (v0 < v) // x0 is bigger than the true value
    {
        do
        {
            x1 *= 0.9;
            v1 = phi1(x1, y);
            //printf("x1=%.16lf, t=%.16lf\n", x1, v1);
        } while (v1 < v);
        *pv = v1;
        return x1;
    }

    // now x0 is smaller than the true value
    do
    {
        x1 *= 1.1;
        v1 = phi1(x1, y);
        //printf("x1=%.16lf, t=%.16lf\n", x1, v1);
    } while (v1 > v);

    *pv = v1;
    return x1;
}
//5. 函数 alpha(v,y) 表示方程 phi1(s,y) - v=0 的解s
double alpha(double v, int y)
{
    assert(y <= 311);

    ALPHA_ST c;
    int i;
    double k, a2;

    y = correct(y);
    c.v = v;
    //printf("v=%.16lf\n", v);

    c.a0 = get_s0_v1(v, y);
    c.v0 = phi1(c.a0, y);

    //printf("a0=%.16lf,v0=%.16f\n", c.a0, c.v0);

    c.a1 = get_x1(y, v, c.a0, c.v0, &(c.v1));
    //printf("a1=%.16lf, v1=%.16lf\n", c.a1, c.v1);

    i = 2;
    while (true)
    {
        k = (c.v1 - c.v0) / (c.a1 - c.a0);
        a2 = c.a1 + (c.v - c.v1) / k;

        c.a0 = c.a1;
        c.v0 = c.v1;
        c.a1 = a2;
        c.v1 = phi1(c.a1, y);
        //printf("\n a[%d]= %.16lf, v[%d]=%.16lf\n", i, c.a1, i, c.v1);
        if (fabs(c.a1 - c.a0) < c.a0 * 1e-9)
            break;
        i++;
    }

    k = (c.v1 - c.v0) / (c.a1 - c.a0);
    a2 = c.a1 + (c.v - c.v1) / k;
    //printf("\n phi1(%.16lf,%d)=%.16lf\n", a2, y, phi1(a2, y));
    return a2;
}

double sm_count(double x, int y)
{
    double s, z, u;
    double r;

    y = correct(y);
    s = alpha(log(x), y);
    z = zeta(s, y);
    u = phi2(s, y);
    r = pow(x, s) * z / (s * sqrt(2*PI * u));
    return r;
}

double sm_count_logx(double logx, int y)
{
    double s, z, u;
    double r;

    y = correct(y);
    s = alpha(logx, y);
    z = zeta(s, y);
    u = phi2(s, y);
    r = s*logx + log(z) - log(s*sqrt(2*PI*u));
    return r;
}

void test()
{
    double x, s, logx;
    int n, y;

    n = 48;
    y = primes[n - 1];
    x = pow(2.0, 64);
    logx = log(x);
    s = alpha(log(x),y);
}

void test2()
{
    double x = pow(2.0, 64.0);
    double r;
    int y;
    for (int i = 4; i <= 48; i += 4)
    {
        y = primes[i - 1];
        r = sm_count(x, y);
        printf("sm_count( 2^64,%d)= %.16lf\n", y,r);
    }
}

void test3()
{
  int y;
  int pow2;
  int n;
  int step;
  for (int j=0;j<2;j++)
  {
     if (j==0)
     {	pow2=64; n=12; step=4;}
     else
     { pow2=4096; n=2; step=4;}
     printf("\n----------------------------------\n");
     for (int i=1;i<=n;i++)
     {
        int y=primes[i*step-1];
   double logx = pow2*log(2.0);
    double logsm= sm_count_logx(logx,y);
    //printf("log(sm_count(2^%d,%d))=%.16lf\n",pow2,y,logsm);
    double r= exp(logsm);
    printf("sm_count(2^%d,%d)=%.16lf\n",pow2,y,r);
     }
  }
}
int main()
{
    test2();
    test3();
    return 0;
}

/* test2 输出结果
sm_count( 2^64,7)= 87145.0247772681177594
sm_count( 2^64,19)= 11286548.3023795727640390
sm_count( 2^64,37)= 211967017.2076353728771210
sm_count( 2^64,53)= 1599684387.4000642299652100
sm_count( 2^64,71)= 7107020116.2066755294799805
sm_count( 2^64,89)= 22998573481.4820442199707031
sm_count( 2^64,107)= 58913795304.6139526367187500
sm_count( 2^64,131)= 129415464204.2234191894531250
sm_count( 2^64,151)= 250108239713.3505859375000000
sm_count( 2^64,173)= 441507555760.2599487304687500
sm_count( 2^64,193)= 726098055915.4711914062500000
sm_count( 2^64,223)= 1128020024926.1152343750000000
*/


/* test3 输出结果
----------------------------------
sm_count(2^64,7)=87145.0247772680886555
sm_count(2^64,19)=11286548.3023795988410711
sm_count(2^64,37)=211967017.2076350152492523
sm_count(2^64,53)=1599684387.4000587463378906
sm_count(2^64,71)=7107020116.2067070007324219
sm_count(2^64,89)=22998573481.4820785522460938
sm_count(2^64,107)=58913795304.6140899658203125
sm_count(2^64,131)=129415464204.2231750488281250
sm_count(2^64,151)=250108239713.3501892089843750
sm_count(2^64,173)=441507555760.2587890625000000
sm_count(2^64,193)=726098055915.4709472656250000
sm_count(2^64,223)=1128020024926.1152343750000000
----------------------------------
sm_count(2^4096,7)=1163396154613.0056152343750000
sm_count(2^4096,19)=884405629135495102464.0000000000000000
*/

 
/* 
n | 真实值 | 计算值  
4; | 85,348 | 87,145.
8; | 11,169,655 | 11,286,548.
12; | 210,499,995 | 211,967,017.
16; | 1,591,363,196 | 1,599,684,387.
20; | 7,077,350,886 | 7,107,020,116.
24; | 22,918,160,232 | 22,998,573,481.
28; | 58,735,924,658 | 58,913,795,304
32; | 129,070,232,683 | 129,415,464,204.
36; | 249,507,947,485 | 250,108,239,713.
40; | 440,539,214,034 | 441,507,555,760.
44; | 724,624,467,775 | 726,098,055,915.
48; | 1,125,880,419,418 | 1,128,020,024,926.
*/