[C语言] 计算光滑数的个数的算法，个人思考，仅供参考

最新推荐文章于 2024-04-17 14:40:43 发布

haitunhaoma

最新推荐文章于 2024-04-17 14:40:43 发布

阅读量317

点赞数

本文链接：https://blog.csdn.net/haitunhaoma/article/details/103741202

版权

本文介绍了一种计算光滑数个数的算法，利用C语言实现。算法涉及了多项数学函数，如zeta函数、phi1和phi2函数，并通过牛顿迭代法求解方程。代码中给出了不同限制条件下的光滑数数量计算示例。

摘要由CSDN通过智能技术生成

#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];