程序師世界是廣大編程愛好者互助、分享、學習的平台,程序師世界有你更精彩!
首頁
編程語言
C語言|JAVA編程
Python編程
網頁編程
ASP編程|PHP編程
JSP編程
數據庫知識
MYSQL數據庫|SqlServer數據庫
Oracle數據庫|DB2數據庫
 程式師世界 >> 編程語言 >> C語言 >> C++ >> C++入門知識 >> poj 3696 The Luckiest number

poj 3696 The Luckiest number

編輯:C++入門知識

這個題很奇葩了。題意是給出個數字L,假如存在一個數K使得L*K = 888...,求888...的最小長度,如果不存在這樣的K,那麼輸出0。
我是什麼思路也沒有了,拖了幾天了,數論搞死我了,只能找答案了。
   我看到個比較靠譜的推法。首先,888...=111...*8=(10^0+10^1+...+10^m-1)*8=(10^m - 1)/9*8,PS:m代表888...的長度。
好吧,終於化成指數了,現在有8*(10^m-1)/9=K*L,最小的m就是我們要求的答案啦。

   方式1:
   => 8 * (10^m-1) = 9 * k * L
   => 8/d*(10^m-1)=9*k*L/d,d=gcd(8,9L)
   => 10^m-1 = 0 % 9 * L / gcd(8, 9L) = 0 % 9*L/gcd(8,L),(由於gcd(8/d,9L/d)=1,那麼10^m-1必然是9*L/d的倍數了)。
   => 10^m = 1 % 9 * L / gcd(8,L)
   方式2:
   => 8*(10^m-1)/9 = 0 % L
   => 8*(10^m-1) = 0 % 9*L(這步的推出,比如x/9 = k*n,那麼x=9*k*n了,顯然成立)
   => 10^m-1 = 0 % 9*L/gcd(9*L,8),假如,d = gcd(9*L,8),那麼有8/d*(10^m-1)=k*9*L/d,因為8/d不可能是9  *L / d
的倍數,所以10^m-1必定是9*L/d的倍數,所以10^m-1 = 0 % 9*L/gcd(9*L,8)),=>,10^m - 1 = 0 % 9 * L / gcd(L, 8),
(因為gcd(9,8)=1)。
   => 10^m = 1 % 9*L/gcd(8,L) 

   至此,2種方式都推出了,10^m = 1 % 9*L/gcd(8,L) 。
   那麼怎麼解答這個問題了,這個就用到了歐拉定理了。令p = 9 * L / gcd(8,L),那麼有10^m = 1 % p。由歐拉定理知,Z*p中所有的
數字a均滿足a^euler(p) = 1 % p。那麼,10只要是p的乘法群中就肯定有解了。如果,10不在Z*p中了,那麼有10^m= 2^m * 5^m。
而且10和p有公告因子2或者5,所以p = 2 * k或者 p = 5 * k,2^m=0%p或者5^m=0%p,那麼10^m就永遠不可能是1%p了。
   綜上所述,要滿足式子a^m=1%p,必須gcd(p,a)=1,即a必須是p的乘法群中的數字。
   現在的問題是求最小的m,由歐拉定理知道a^euler(p)=1%p,m再大就開始循環了。但是m可能會更小。比如,我們現在知道最小的m
是min,那麼有a^min=1%p,因為要滿足a^euler(p)=1%p,那麼a^euler(p)肯定能變換成(a^min)^k,至於k是多少就不知道了,當然
也可以求出來。那麼min就是euler(p)的一個因子,而且是最小的一個滿足a^min=1%p的因子了。
   現在就可以通過枚舉euler(p)的因子,找到最小的因子min滿足式子a^min = 1 % p就能解決本問題了。
   注意求a^m%p肯定是通過算法導論上面那種方法的,O(32)或者O(64)的復雜度,還有a*b%m也需要自己模擬,因為可能a*b就溢出了。
   代碼如下,貌似代碼還可以通過其它的改進加快速度。

#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <string.h>
using namespace std;
typedef long long INT;

//10^m = 1 % (9*L / gcd(8, L)),求最小m
//p = 9 * L / gcd(8,L)
//gcd(p,10) != 1則p有2或者5的因子,2^m=1%p或者
//5^m=1%p無解,原式無解
//if(p)素數,m=euler(p) = p - 1
//否則,m一定是euler(p)的最小滿足等式的因子
//因為(10^m)^n = 10^euler(p) = 1%p
INT gcd(INT a, INT b)
{
    if (a < b)swap(a, b);
    while (b)
    {
        INT t = a;
        a = b;
        b = t % b;
    }
    return a;
}

INT Euler(INT nN)
{
    INT nAns = 1;
    INT nMax = sqrt((double)nN) + 1;
    for (INT i = 2; i <= nMax; ++i)
    {
        if (nN % i == 0)
        {
            nAns *= i - 1;
            nN /= i;
            while (nN % i == 0)
            {
                nAns *= i;
                nN /= i;
            }
        }
    }
    if (nN != 1)nAns *= nN - 1;
    return nAns;
}

INT MultMod(INT a, INT b, INT mod)
{
    INT ans = 0;
    while (b)
    {
        if (b & 1)
        {
            ans = (ans + a) % mod;
        }
        a = (2 * a) % mod;
        b >>= 1;
    }
    return ans;
}

INT ExpMod(INT base, INT exp, INT mod)
{
    INT ans = 1;
    base %= mod;
    while (exp)
    {
        if (exp & 1)
        {
            ans = MultMod(ans, base, mod);
        }
        base = MultMod(base, base, mod);
        exp >>= 1;
    }
    return ans % mod;
}

INT GetAns(INT p)
{
    INT u = Euler(p);
    INT nMax = sqrt((double)u) + 1;
    INT nAns = u;
    for (INT i = 1; i <= nMax; ++i)
    {
        if (u % i == 0)
        {
            if (ExpMod(10, i, p) == 1)
            {
                nAns = i;
                break;
            }
            if (ExpMod(10, u / i, p) == 1)
            {
                nAns = min(nAns, u / i);
            }
        }
    }
    return nAns;
}

int main()
{
    INT nL;
    INT nCase = 1;
   
    while (scanf("%I64d", &nL), nL)
    {
        INT p = 9 * nL / gcd(nL, 8);
        if (gcd(p, 10) != 1)
        {
            printf("Case %I64d: 0\n", nCase++);
            continue;
        }
        printf("Case %I64d: %I64d\n", nCase++, GetAns(p));
    }
   
    return 0;


  1. 上一頁:
  2. 下一頁:
Copyright © 程式師世界 All Rights Reserved