#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <limits.h>
#include <float.h>
#define _CRT_SECURE_NO_WARNINGS
#define LN10 2.302585092994045684018
#define RLN10 0.43429448190325182765
#define MAX_LEN 128
// 1+2 +3 + ... q = (q+1)q/2
// 1^2+2^2 +3^2 + ... q^2 = q(q+1)(2q+1)/6
// 1^3+2^3 +3^3 + ... q^3 = (q(q+1)/2)^2
typedef struct
{
double v; // p/q
double e1; // 1次项 (1+2+3 + ...q)/LN10= (q+1)*q/(2*LN10)
double e2; // 2次项 (1^2+2^2+3^2 + ...q^2)/(2*LN10)=q*(q+1)*(2q+1)/(12*LN10)
double e3; // 3次项 (1^3+2^3+3^3 + ...q^3)/(3*LN10)= ((q+1)q)^2/(12LN10)
int32_t p;
int32_t q;
int idx;
} FRAC;
typedef struct _dblArray
{
double lgn; // log10(n)- floor( log10(n) )
int n;
FRAC ls[MAX_LEN];
FRAC hs[MAX_LEN];
}DBL_ARRAY;
inline double get_double(FRAC *pf)
{
return (double)pf->p / (double)pf->q;
}
// 1+2 +3 + ... q = (q+1)q/2
// 1^2+2^2 +3^2 + ... q^2 = q(q+1)(2q+1)/6
// 1^3+2^3 +3^3 + ... q^3 = (q(q+1)/2)^2
inline void update_FRAC get_double(FRAC *p)
{
//double e1; // 1次项 (1+2+3 + ...q)/LN10= (q+1)*q/(2*LN10)
//double e2; // 2次项 (1^2+2^2+3^2 + ...q^2)/(2*LN10)=q*(q+1)*(2q+1)/(12*LN10)
//double e3; // 3次项 (1^3+2^3+3^3 + ...q^3)/(3*LN10)= ((q+1)q)^2/(12LN10)
double f = p->q;
double e1 = f * (f + 1.0) * 0.5;
p->e2 = f * (f + 1.0) * (f * 2 + 1) / (12.0 * LN10);
p->e3 = (e1 * e1) / (3.0 * LN10);
p->e1 = e1 / LN10;
p->v= (double)(p->p) / (double)(p->q);
}
void create_DblArray ( DBL_ARRAY *dblArray, int n, double lgn)
{
FRAC *pf1;
FRAC *pf2;
dblArray->lgn = lgn;
dblArray->n = n;
pf1 = dblArray->ls + 0;
pf2 = dblArray->hs + 0;
pf1->p = 1;
pf1->q = (int32_t)ceil(1.0 / lgn);
pf1->idx = 0;
update_FRAC(pf1);
pf2->p = 1;
pf2->q = (int32_t)floor(1.0 / lgn);
pf2->idx = 0;
update_FRAC(pf2);
}
// (lg(n+d)!)的小数部分
// (lg(n + d)!)
// ~= lg(n) + d*lg(n) + d*(d+1)/(2n*ln(10)
inline double get_lg_new_n(double lgn, FRAC fr, double rk)
{
double det = lgn - get_double(fr); //
double r = det * (double)(fr.q) + (double)(fr.q * (fr.q + 1)) * rk;
return r - round(r);
}
//首次算出欠周期和过周期
//当Q很大时或者欠周期>0时，使用常规方法求得步进
// 在计算序列时，一旦不收敛时（欠周期时得到的值>0或者过周期得到的值不再减少），直接退出
// 当欠时，使用过周期，取最后一个累积和大于0的数
// 当过时，使用欠周期，取最后一个累积和小于0的数
void farey_search(uint64_t n, double f_lgn, double det)
{
FRAC low, high, mid;
uint32_t best_q;
double rk = 1.0 / ((double)n * 2.0 * LN10);
double big_lg_nn, small_lg_nn; //big lg(new_n) and small lg(new_n)
double last_det, new_det, lg_nn, best_diff;
bool isfirst;
int s_det, s_ndet;
low.p = 1;
low.q = ceil(1.0 / f_lgn);
high.p = 1;
high.q = floor(1.0 / f_lgn);
if (low.p >= 65535 || get_lg_new_n(f_lgn, high, rk) >= 0)
{
printf("Need use common method to do this when n=%lu", n);
return;
}
printf("original det =%.16lf\n", det);
big_lg_nn = get_lg_new_n(f_lgn, low, rk); //过周期
small_lg_nn = get_lg_new_n(f_lgn, high, rk); //欠周期
if (det < 0) //使用过周期
{
best_diff = det + big_lg_nn;
best_q = low.q;
}
else
{
best_diff = det + small_lg_nn;
best_q = high.q;
}
isfirst = true;
last_det = det;
while (true)
{
mid.p = low.p + high.p;
mid.q = low.q + high.q;
printf("mid.q=%d\n", mid.q);
lg_nn = get_lg_new_n(f_lgn, mid, rk);
new_det = lg_nn + det;
if (get_double(mid) >= f_lgn) //是上界，欠周期
{
if (lg_nn >= 0.0) // 欠周期时得到的值 >= 0
{
printf("Break from %s:%d\n", __FILE__, __LINE__);
break;
}
if (det > 0)
{
if (isfirst || fabs(new_det) < fabs(last_det))
{
printf("run to %s:%d\n", __FILE__, __LINE__);
best_q = mid.p;
best_diff = new_det;
}
else
{
printf("Break from %s:%d\n", __FILE__, __LINE__);
break;
}
}
high = mid;
small_lg_nn = lg_nn;
}
else //是下界，过周期
{
if (lg_nn >= big_lg_nn) //过周期得到的值不再减少
{
printf("break from %s:%d\n", __FILE__, __LINE__);
break;
}
if (det < 0)
{
if (isfirst || fabs(new_det) < fabs(last_det))
{
printf("run to %s:%d\n", __FILE__, __LINE__);
best_q = mid.p;
best_diff = new_det;
}
else
{
printf("break from %s:%d\n", __FILE__, __LINE__);
break;
}
}
low = mid;
big_lg_nn = lg_nn;
}
if (mid.q >= 65535)
{
printf("run to %s:%d\n", __FILE__, __LINE__);
break;
}
last_det = new_det;
isfirst = false;
}
printf("best q is %u,new det =%.16lf\n", best_q, best_diff);
}
int main(int argc, char *argv[])
{
// lgpi = log(3.1415926535898) / log(10)
// n = 1544607046599
// lgn = log(n) / log(10)
// lg_fn = lngamma(n + 1) / log(10)
// f_lg_fn = lg_fn - floor(lg_fn)
// f_lgn = lgn - floor(lgn)
// det = lgpi - f_lg_fn
// 0.188818011786880272870
// - 2.0713896587241805909e-14
// 0.1888180120936348488209645817901717275944058254028858993817802728002129534568670651162967209435905600
// -0.0004510269743969044278283267303117723666536375937437949826689256603143119413249000389634605276822730510
//uint64_t n = 1544607046599;
//double f_lgn = 0.188818011786880272870;
//double det = -2.0713896587241805909e-14;
uint64_t n = 1544607046599 + 1091;
double f_lgn = 0.188818012093634848820964581790;
double det = -0.000451026974396904427828326734;
printf("n=?");
scanf("%lu", &n);
printf("f_lgn=?");
scanf("%lf", &f_lgn);
printf("det=?");
scanf("%lf", &det);
farey_search(n, f_lgn, det);
return 0;
}