[HDOJ3037 Saving Beans]Lucas定理

【题目大意】
求C(n+m,n) % p的值。
保证p是素数。
【算法分析】
师傅说存在一个叫Lucas定理。
他在维基百科上:http://en.wikipedia.org/wiki/Lucas%27_theorem
除了这个Lucas定理以外,还要用到快速幂和另外一个不知名定理:
若b与p互质,则(a/b) = a * b^(p-2) (mod p)。
【其它】算阶乘的时候弄错几次,要开int64,不然会乘到溢出的。
【CODE】
#include #include #include #include using namespace std;
typedef __int64 lld;
lld p;
lld jc[105555];

void init(){
jc[0]=1;
for (lld i=1;i<=p;i++) jc[i]=jc[i-1]*i%p;
}

lld Pow_mod(lld x,lld cf){
if (cf==0) return 1;
lld res=Pow_mod(x,cf/2);
res=res*res%p;
if (cf%2) res=res*x%p;
return res;
}

lld C(lld m,lld n){
if (n>m) return 0;
return jc[m]*Pow_mod(jc[n]*jc[m-n]%p,p-2)%p;
}

lld solve(lld m,lld n){
if (n==0) return 1;
return C(m%p,n%p)*solve(m/p,n/p)%p;
}

int main(){
lld Tc,i,m,n;
cin >> Tc;
for (i=0;i cin >> n >> m >> p;
init();
cout << solve(m+n,n) << "n";
}
return 0;
}

[FOJ1759 Super A^B mod C]某数论公式的应用

【题目】
求A^B % C
A,C<=1e9。
B<=10^1000000
【算法分析】
师傅说有一个公式
A^x=A^(x % phi(c) + phi(c) )    (mod c)
然后直接套。
【其它】1A
【CODE】
#include #include #include #include using namespace std;
typedef long long lld;
const lld N=1005555;
lld A,B,C,phi;
char str[N];

void makephi(){
lld x=C,i;
phi=C;
for (i=2;i*i<=x;i++)
if (x%i==0){
phi=phi/i*(i-1);
while (x%i==0) x/=i;
}
if (x!=1) phi=phi/x*(x-1);
}

void String_to_Number(){
bool flag=false;
lld i,j,Len=strlen(str+1);
B=0;
for (i=1;i<=Len;i++){
B=B*10+str[i]-‘0’;
if (B>=phi){
flag=true;
B%=phi;
}
}
if (flag) B+=phi;
}

lld Pow_mod(lld x,lld cf){
if (!cf) return 1;
lld res=Pow_mod(x,cf/2);
res*=res; res%=C;
if (cf%2) res*=x;
res%=C;
return res;
}

int main(){
while (scanf("%I64d %s %I64d",&A,str+1,&C)!=EOF){
makephi();
String_to_Number();
printf("%I64dn",Pow_mod(A,B));
}
return 0;
}

[SPOJ LCMSUM]数论

【题目大意】
求LCM(1,N)+LCM(2,N)+LCM(3,N)+…+LCM(N,N)
【算法分析】
ans=n*(1/gcd(1,n)+2/gcd(2,n)+…+n/gcd(n,n))
其中显然gcd(x,n)都是n的约数。
然后我们都知道n的约数的个数是n^0.5级别的。
于是可以考虑把分母相同的合并。
然后设x为n的一个约数,那么(kx,n)=x 当且仅当满足 (k,n/x)=1
然后再有不等式kx<=n 即 k<=n/x。
然后。。。
k<=n/x   (k,n/x)=1     想到了什么?
欧拉函数~
设t=n/x
但是现在要求的是所有满足  (k<=t)&& (k,t)=1  的k的和。
然后我们来证一个性质。
那个和=t*phi(t)/2。

假设:gcd(x,t)=g && g>1
那么gcd(t-x,t)%g=0
为什么?
∵t%g=0   x%g=0  ∴ (t-x) % g=0。
所以可以利用反证法得出,假如(x,t)=1那么 (t-x,t)=1
于是互质的数都是成对的,且他们的和是t。
所以总体来说可以分成phi(t)/2组,每组的和为t,所以总和为t*phi(t)/2。

但是如果t为偶数的话,有可能出现t/2算两次的情况之类。
但是显然,t为偶数时,t/2和t不互质,不会算进欧拉函数里,所以不用考虑。

然后最后的做法是:
1、
用筛法把1e6的素数和欧拉函数都筛出来。 这一步的时间复杂度是O(1e6 lg 1e6)。
这个时间复杂度可能很多人会说= =,但是至少我在网上看到的筛法好像都是O(n lg n)的复杂度的。
虽然他们叫线性筛,我也一直认为是线性的,直到APIO听了CDQ讲课,其中有一个典型的复杂度分析就是我这种写法的筛。最后的结论是O(n lg n)。。。
2、
先求n的所有约数,然后利用上面的性质求和。
但是这题异常猥琐,时间非常紧,要先用筛出的素数对n分解质因数,然后再DFS求约数。
一开始我过于信任SPOJ的机子,直接一个for过去求约数,复杂度O(n^0.5)。
然后除了极限数据,我自己的机子卡着5s过。但是上面就不行了。

用了分解质因数以后,就由于那些素数的小步伐跳跃,卡过了。。。~

【时间复杂度】
O(1e6 lg 1e6+ ( T(n^0.5) + lg n )*TestCaseNum)。
其中T(x)函数表示<=x的数里的素数个数。
然后第二步如果是我本来那种直接暴力找约数的话,那么复杂度就是
O(1e6 lg 1e6+ ( n^0.5 + lg n)*TestCaseNum)。
最终,由于TestCaseNum把常数无限放大,我的暴力悲剧地TLE了。。。
【其它】
在HWG的博客看到这题。。。于是做了一下,没想到被常数卡成SB了。。。
Orz,一开始以为A|B是AB互质的意思,现在已经修正。。。
【CODE】
#include #include #include typedef long long lld;
const int N=1000001;
int n,Tc,pt,Lt,sum;
bool prime[N];
int phi[N],pl[N],L[N],Num[N];
lld ans;
char ss[200];

void init(){
memset(prime,true,sizeof(prime));
pt=0;
int i,j;
for (i=1;i for (i=2;i for (pl[++pt]=i,j=i;j phi[j]=phi[j]/i*(i-1);
prime[j]=false;
}
}

void Div(int x){
Lt=0;
for (int i=1;(pl[i]*pl[i]<=n) && (x!=1);i++){
if (x%pl[i]==0) {L[++Lt]=pl[i]; Num[Lt]=0;} else continue;
while (x%pl[i]==0){
Num[Lt]++;
x/=pl[i];
}
}
}

void dfs(int dep){
if (dep>Lt){
ans+=((lld)(phi[sum])*sum)>>1;
if (sum*sum!=n) ans+=((lld)(phi[n/sum])*(n/sum))>>1;
return;
}
int tmp=sum;
dfs(dep+1);
for (int i=1;i<=Num[dep];i++){
sum*=L[dep];
if (sum>n/sum) break;
dfs(dep+1);
}
sum=tmp;
}

void Write(){
int tt=0;
while (ans){
ss[tt++]=ans%10+’0′;
ans/=10;
}
for (int i=tt-1;i>=0;i–) putchar(ss[i]);
putchar(‘n’);
}

void solve(){
phi[1]=2;
sum=1;
ans=0;
Div(n);
dfs(1);
ans*=n;
printf("%lldn",ans);
}

void Read(int &x){
char ch;
for (ch=getchar();!(ch>=’0′ && ch<='9');ch=getchar());
x=0;
for (;ch>=’0′ && ch<='9';ch=getchar())
x=x*10+ch-‘0’;
}

int main(){
init();
Read(Tc);
int i;
for (i=1;i<=Tc;i++){
Read(n);
solve();
}
return 0;
}

[SGU 217 && HDOJ 3310]微积分

【题目大意】
给两个圆柱体,他们的中轴线共面且垂直。而且这两个圆柱体的高是无限大的。
然后分别给出他们底面的半径r1和r2。
求他们相交的体积。
【算法分析】
按X或Y轴积分怎么搞都搞不过。。。。我倒。。。
按Z轴一积就过了。
还是两个都啰嗦一下吧。
按X或Y轴切得话,切面就是这种形状:“|○|” 然后分类讨论一下那两个竖线夹得地方算面积就可以了。
按Z轴切得话,切面就是”井“字形。然后算矩形面积即可。
【其它】
估计是按X或Y轴切时算面积涉及到arccos和π,所以使误差加大了。。。
【CODE】
#include #include #include #include #define sqr(x) ((x)*(x))
const double pi=acos(-1);
double r1,r2,ans,len,eps,tmp;

inline double f(double x){return sqrt((sqr(r1)-sqr(x))*(sqr(r2)-sqr(x)));}

void work(){
ans=0;
for (double x=0;x+eps ans+=(f(x)+f(x+eps)+f(x+eps/2))/3*eps;
ans*=8;
printf("%.2lfn",ans);
}

int main(){
int Tc;
scanf("%d",&Tc);
for (int i=1;i<=Tc;i++){
scanf("%lf%lf",&r1,&r2);
if (r1>r2){tmp=r1; r1=r2; r2=tmp;}
eps=r1/1500000;
work();
}
}

[ZSOI 三核苷酸]求一个特殊的方差

【题目】

Description

  三核苷酸是组成DNA序列的基本片段。具体来说,核苷酸一共有4种,分别用’A’,’G’,’C’,’T’来表示。而三核苷酸就是由3个核苷酸排列而 成的DNA片段。三核苷酸一共有64种,分别是’AAA’,’AAG’,…,’GGG’。给定一个长度为L的DNA序列,一共可以分辨出(L-2)个三核 苷酸。现在我们想用一些统计学的方法来进行一些分析,步骤如下:

给定一个DNA序列,请你计算出它的方差。

Input

  输入包含多组测试数据。第一行包含一个正整数T,表示测试数据数目。每组数据包含一个由’A’,’G’,’C’,’T’组成的字符串,代表要统计的 DNA序列。DNA序列的长度大于等于3且不会超过100000。

Output

  对每组测试数据,输出一行答案,为一个保留6位精度的实数,代表S2的值。如果你的答案和标准答案的“相对误差”小于1e-8,你的答案会被视为正确 的答案。

Sample Input

1
ATATATA

Sample Output

0.750000

【算法分析】
由于只有64种那个XX,所以开一个大小为64的桶,并将之与4进制数对应起来。
然后就是类似dp那样,根据前面的信息,推出新的信息。
【其他】
由于精度过紧,而我又不会用C++的long double类型的scanf读入和printf输出。。。
于是果断Pascal。。。
【CODE】
const
inf=’tri.in’;
ouf=’tri.out’;

var
Tc,n:longint;
p,total:array [0..64] of extended;
num,ss:array [0..64] of extended;
str:ansistring;
a:array [1..200000] of longint;

procedure init;
var
i:longint;
begin
for i:=1 to n do
case str[i] of
‘A’:a[i]:=0;
‘T’:a[i]:=1;
‘G’:a[i]:=2;
‘C’:a[i]:=3;
end;
end;

procedure solve;
var
i,now:longint;
sum,times,pj:extended;
begin
sum:=0; times:=0;
fillchar(num,sizeof(num),0);
fillchar(total,sizeof(total),0);
fillchar(p,sizeof(p),0);
for i:=n-2 downto 1 do
begin
now:=a[i]*4+a[i+1];
now:=now*4+a[i+2];
if p[now]=0 then
begin
total[now]:=1;
p[now]:=i;
num[now]:=0;
end
else
begin
sum:=sum+total[now]*(p[now]-i)+num[now];
times:=times+total[now];

num[now]:=total[now]*(p[now]-i)+num[now];
total[now]:=total[now]+1;
p[now]:=i;
end;
end;
if times=0 then
begin
sum:=0;
writeln(sum:0:6);
exit;
end;
pj:=sum/times;
fillchar(num,sizeof(num),0);
fillchar(total,sizeof(total),0);
fillchar(p,sizeof(p),0);
fillchar(ss,sizeof(ss),0);
sum:=0;
for i:=n-2 downto 1 do
begin
now:=a[i]*4+a[i+1];
now:=now*4+a[i+2];
if (p[now]=0) then
begin
total[now]:=1;
p[now]:=i;
num[now]:=0;
ss[now]:=0;
end
else
begin
num[now]:=num[now]
+sqr(p[now]-i)*(total[now]-1)
+2*ss[now]*(p[now]-i)+sqr(p[now]-i-pj);

ss[now]:=ss[now]
+(total[now]-1)*(p[now]-i)
+p[now]-i-pj;

total[now]:=total[now]+1;

p[now]:=i;

sum:=sum+num[now];
end;
end;
writeln(sum/times:0:6);
end;

procedure main;
var
i:longint;
begin
readln(Tc);
for i:=1 to Tc do
begin
readln(str);
n:=length(str);
init;
solve;
end;
end;

begin
assign(input,inf);
assign(output,ouf);
reset(input);
rewrite(output);
main;
close(input);
close(output);
end.