#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <stdint.h>
#include <math.h>
#include <mpfr.h>
#include <mpreal.h>
#include <time.h>
#include <inttypes.h>
#define LG2 0.30102999566398119521
#define RLN2 1.44269504088896340735
#define LOG_PIx2 1.83787706640934548356
time_t g_last_time;
typedef mpfr::mpreal real;
void printTab(int c)
{
if (c>1024)
c=1024;
printf("{\n ");
for (int i=0;i<c;i++)
{
int k= int(floor(LG2*i));
if ( k==0)
printf("1");
else
printf("1e%d",-k);
if ( i<c-1)
printf(",");
if ( i %8==7 || i==c-1)
{
printf("\n");
if ( i != c-1)
printf(" ");
}
}
printf("};\n");
}
typedef union
{
double f;
uint32_t i[2];
}Pack;
#define CHEKC_MIN_MAX(a,b,min,max) { if (a < min) min=a; if (b < min) min=b; if (a > max) max=a; if (b > max) max=b; }
double pow2log10Tab[]=
{
1,1,1,1,1e-1,1e-1,1e-1,1e-2,
1e-2,1e-2,1e-3,1e-3,1e-3,1e-3,1e-4,1e-4,
1e-4,1e-5,1e-5,1e-5,1e-6,1e-6,1e-6,1e-6,
1e-7,1e-7,1e-7,1e-8,1e-8,1e-8,1e-9,1e-9,
1e-9,1e-9,1e-10,1e-10,1e-10,1e-11,1e-11,1e-11,
1e-12,1e-12,1e-12,1e-12,1e-13,1e-13,1e-13,1e-14,
1e-14,1e-14,1e-15,1e-15,1e-15,1e-15,1e-16,1e-16,
1e-16,1e-17,1e-17,1e-17,1e-18,1e-18,1e-18,1e-18,
1e-19,1e-19,1e-19,1e-20,1e-20,1e-20,1e-21,1e-21,
1e-21,1e-21,1e-22,1e-22,1e-22,1e-23,1e-23,1e-23,
1e-24,1e-24,1e-24,1e-24,1e-25,1e-25,1e-25,1e-26,
1e-26,1e-26,1e-27,1e-27,1e-27,1e-27,1e-28,1e-28,
1e-28,1e-29,1e-29,1e-29,1e-30,1e-30,1e-30,1e-31,
1e-31,1e-31,1e-31,1e-32,1e-32,1e-32,1e-33,1e-33,
1e-33,1e-34,1e-34,1e-34,1e-34,1e-35,1e-35,1e-35,
1e-36,1e-36,1e-36,1e-37,1e-37,1e-37,1e-37,1e-38
};
int min_prec( uint64_t n, int digLen)
{
double logfac;
logfac = 0.5*LOG_PIx2 + (n+0.5)*log((double)n)-n;
logfac *= RLN2;
//printf("log2(%" PRIu64 "!)=%.16lf\n", n,logfac);
double minPrec= log(logfac)/log(2)+ digLen*3.33;
return (int)minPrec;
}
double get_fac( uint64_t n)
{
double r=1.0;
if ( n<171)
{
for (uint32_t i=1;i <= n;i++)
r*= (double)i;
uint32_t *p= (uint32_t*)(&r);
int e2 = (p[1]>> 20) - 1023;
int e10= int(LG2*e2)-1;
r /= pow(10.0,(e10>0?e10:0) );
while ( r >= 10.0)
r*=0.1;
return r;
}
mpfr::mpreal::set_default_prec(min_prec(n,20));
real ln10= mpfr::log(10);
r= (double)mpfr::pow(10,mpfr::frac(mpfr::lngamma(n+1)/ln10) );
return r;
}
bool check_fac( uint64_t n, uint64_t pi, int digLen)
{
char fmt[32];
real lowBound,hiBound,r;
mpfr::mpreal::set_default_prec( min_prec(n,digLen+8) );
lowBound = pi; lowBound /= mpfr::pow(10.0, digLen-1);
hiBound = pi+1; hiBound /= mpfr::pow(10.0, digLen-1);
real ln10 =mpfr::log(10);
r= (double)mpfr::pow(10,mpfr::frac(mpfr::lngamma(n+1)/ln10) );
if (r >= lowBound && r< hiBound)
{
uint64_t log10n = (uint64_t)(mpfr::rint_floor(mpfr::lngamma(n+1)/ln10)) ;
printf("%" PRIu64 "! =", n);
sprintf(fmt,"%%.%dR*f", digLen+9 );
mpfr_printf(fmt, MPFR_RNDN, r);
printf("* 10^%" PRIu64 "\n", log10n);
return true;
}
else
return false;
}
bool search_fac(uint64_t base, uint64_t count, uint64_t pi, int digLen, uint64_t &lastBase )
{
Pack p[2];
double lowBound,hiBound,det;
double min,max;
int k0,k1;
det= DBL_EPSILON * 4.0 * count;
lowBound = pi / pow(10.0, digLen-1); lowBound -= lowBound *det;
hiBound = (pi+1) / pow(10,digLen-1); hiBound += hiBound * det;
//printf("lowBound= %.16lf, hiBound=%.16lf\n",lowBound,hiBound);
//min=100; max=0.01;
p[0].f=p[1].f= get_fac(base-1);
uint64_t end= base+count;
//printf("search %" PRIu64 " to %" PRIu64 ",pi=%" PRIu64 ",digLen=%d\n", base,end,pi,digLen);
uint64_t i=base;
double fi=(double)base;
for (; i < end; i +=2, fi += 2.0)
{
p[0].f = p[1].f * fi;
p[1].f = p[0].f * (fi+1.0);
k0 = (p[0].i[1]>> 20) - 1023;
k1 = (p[1].i[1]>> 20) - 1023;
p[0].f *= pow2log10Tab[k0];
p[1].f *= pow2log10Tab[k1];
//CHEKC_MIN_MAX( p[0].f, p[1].f, min,max)
if ( p[0].f < hiBound && p[0].f >= lowBound)
{
if ( check_fac( i, pi, digLen) )
{ lastBase= i; return true; }
}
if ( p[1].f < hiBound && p[1].f >= lowBound)
{
if ( check_fac( i+1, pi, digLen) )
{ lastBase= i+1; return true; }
}
}
//printf("min=%lf, max=%lf\n", min,max);
time_t now; time(&now);
if ( difftime(now,g_last_time) >= 60) // 每分钟显示一次
{
g_last_time=now;
struct tm* p = localtime(&now);
char timeDir[64];
sprintf(timeDir, "%04d-%02d-%02d %02d:%02d:%02d",
p->tm_year+1900, p->tm_mon+1, p->tm_mday,
p->tm_hour, p->tm_min, p->tm_sec);
printf("%s : up to %" PRIu64 "\n", timeDir, end-1);
}
return false;
}
#define MM_TIMES 10000000 // 最大接受1000万次累乘
void search_all()
{
struct pair { uint64_t n; int l; };
pair arr[] =
{
{ 3,1 },
{ 31,2 },
{ 314,3 },
{ 3141, 4 },
{ 31415, 5 },
{ 314159, 6 },
{ 3141592, 7 },
{ 31415926, 8 },
{ 314159265, 9 },
{ 3141592653, 10 },
{ 31415926535, 11 },
{ 314159265358, 12 },
{ 3141592653589, 13 },
{ 31415926535897, 14 },
{ 314159265358979, 15 },
{ 3141592653589793, 16 },
{ 31415926535897932, 17 },
{ 314159265358979323, 18 },
{ 3141592653589793238, 19 }
};
time(&g_last_time);
uint64_t base=1,lastBase=1;
for (int i=0; i< sizeof(arr)/sizeof(arr[0]); i++ )
{
//printf("lastbase is %" PRIu64 "\n",lastBase);
base=lastBase;
while (true)
{
if ( search_fac(base,MM_TIMES,arr[i].n,arr[i].l,lastBase) )
break;
base+=MM_TIMES;
}
}
}
int main(int argc, char* argv[])
{
search_all();
return 0;
}