exp(pi * sqrt(n))が整数に近くなるnを探す
Posted feedbacks - C#
doubleだと桁が足りなかったのでdecimalにした。そのせいでExpやSqrtを自分で書かないといけなかった。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 | static void Main(string[] args)
{
for (int n = 1; n < 200; n++)
{
decimal f = Func(n);
if (IsAlmostInteger(f, 4))
Console.WriteLine("{0} {1:f6}", n, f);
}
Console.ReadKey();
}
static decimal Func(int i)
{
return DecimalExpQuick(3.141592653589793238462643383279M * DecimalSqrt((decimal)i));
}
static decimal DecimalExpQuick(decimal d)
{
int n = 0;
decimal p = 1.0M;
decimal t = d;
do
{
n++;
p *= 2;
t /= 2;
}
while (t >= 0.1M);
decimal et = DecimalExp(t);
decimal ex = 1;
for(int i = 0; i < p; i++)
ex *= et;
return ex;
}
static decimal DecimalExp(decimal t)
{
int i = 0;
decimal et = 0.0M;
decimal ai = 1.0M;
do
{
et += ai;
i++;
ai = ai * (t / i);
}
while (!NearlyEquals(ai, 0.0M));
return et;
}
static decimal DecimalSqrt(decimal d)
{
decimal x = 1.0M;
decimal betterX = 0.0M;
while(true)
{
betterX = SqrtBailey(x, d);
if (NearlyEquals(x, betterX))
return betterX;
x = betterX;
}
}
static decimal SqrtBailey(decimal x, decimal d)
{
decimal fx = x * x - d;
decimal dfx = 2.0M * x;
return x - (fx / (dfx - (fx / dfx)));
}
static bool NearlyEquals(decimal d1, decimal d2)
{
int dights = IntDight(d1);
decimal delta = 1e-28M * (decimal)Math.Pow(10, dights);
return ((d1 - d2 < delta) && (d2 - d1 < delta));
}
static int IntDight(decimal d)
{
if (d < decimal.Zero)
d = decimal.Negate(d);
return IntAbsDight(d);
}
static int IntAbsDight(decimal d)
{
int i = 0;
while (d >= decimal.One)
{
d /= 10.0M;
i++;
}
return i;
}
static bool IsAlmostInteger(decimal d, int dights)
{
decimal p = (decimal)Math.Pow(10, dights);
decimal d1 = decimal.Round(d) * p;
decimal d2 = decimal.Round(d * p);
return ((d1 - d2 <= 1.0M) && (d2 - d1 <= 1.0M));
}
|


herumi
#3416()
Rating0/2=0.00
Pythonで34行のスクリプトを書いて得られた出力の例が下のようになります。
この問題は光成さんに教えて頂いた e^{π*sqrt{163}}≒26253741640768744 が元になっています。ご協力ありがとうございました。[ reply ]