bzoj 2468 jzoj 2197. 【中山市选2010】三核苷酸

Description

三核苷酸是组成DNA序列的基本片段。具体来说，核苷酸一共有4种，分别用’A’，’G’，’C’，’T’来表示。而三核苷酸就是由3个核苷酸排列而成的DNA片段。三核苷酸一共有64种，分别是’AAA’，’AAG’，…，’GGG’。给定一个长度为L的DNA序列，一共可以分辨出（L-2）个三核苷酸。现在我们想用一些统计学的方法来进行一些分析，步骤如下：
1.对于这（L-2）个三核苷酸，我们从左到右给予编号，分别为1到L-2。
2.从这（L-2）个三核苷酸挑选一对出来，一共有(L-2)*(L-3)/2种可能。如果某一对三核苷酸是一样的，我们就记录他们之间的距离。他们之间的距离定义为他们的编号之差。
3.根据我们所记录的“样本数据”，我们现在需要计算样本数据的方差。方差的计算公式是S2=[(x1-X) 2+(x2-X) 2+…+(xn-X)2]/n, X=(x1+x2+…+xn)/n。如果样本的大小n=0，那么我们认为S2=X=0。

例如，我们要统计DNA序列’ATATATA’：
1. 为三核苷酸编号. L1: ATA, L2:TAT, L3:ATA, L4:TAT, L5:ATA.
2. (L1,L3)=2, (L1,L5)=4, (L3,L5)=2, (L2,L4)=2. 所以样本数据是2,4,2,2.
3. 样本数据平均值X=(2+4+2+2)/4=2.5.
方差S2=[(2-2.5)2+(4-2.5) 2+(2-2.5)2+(2-2.5)2]/4=0.75.
给定一个DNA序列，请你计算出它的方差。

Input

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

Output

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

Sample Input

1
ATATATA

Sample Output

0.750000

分析：设ave为平均数,k为有多少对,t[i]为有每对的值对,答案是

ANS=(ave*ave*k+sum(t[i]^2)-2*ave*sum(t[i]))/k

但是我们不用得到每个t[i],显然只要知道所有t[i]的和以及平方和。

设s[i]为每对的位置,f[i]为3-i，有多少个片段与i的相同，

t[i]=s[i]-s[j] {j为前面和该片段相同的位置}

显然可以记录前面的片段和一次计算，c为上一个位置

t[i]=s[i]*f[c]-sum(s[j])

至于平方和，可以考虑

t[i]^2=(s[i]-s[j])^2

然后套路一波。

注意:
平方和会炸longlong，一定要把k先算出来。
可能一个对都没有，(例如：GGG)，此时应该输出0.000000

代码：

const
 maxn=100001;
var
 t,n:longint;
 f,sums,last:array [0..maxn] of int64;
 square_sums:array [0..maxn] of extended;
 ls:array [0..63] of longint;
 x,y,z,c:longint;
 sumt,sum:int64;
 ave,ans,square_sumt:extended;
 s:ansistring;

function getnum(ch:char):longint;
 begin
  if ch='A' then exit(0);
  if ch='G' then exit(1);
  if ch='C' then exit(2);
  if ch='T' then exit(3);
 end;

procedure main;
 var i,j,k:longint;
begin
 readln(s);
 n:=length(s);
 fillchar(sums,sizeof(sums),0);
 fillchar(square_sums,sizeof(square_sums),0);
 fillchar(ls,sizeof(ls),0);
 sumt:=0; square_sumt:=0; sum:=0;
 x:=getnum(s[1]);
 y:=getnum(s[2]);
 for i:=3 to n do
  begin
   z:=getnum(s[i]);
   c:=ls[x*16+y*4+z];
   f[i]:=f[c]+1;
   sums[i]:=sums[c]+i;
   if f[c]>0 then
     sumt:=sumt+i*f[c]-sums[c];
   last[i]:=c;
   ls[x*16+y*4+z]:=i;
   x:=y; y:=z;
  end;
 for j:=0 to 63 do
  sum:=sum+f[ls[j]]*(f[ls[j]]-1) div 2;
 if sum=0 then begin writeln('0.000000'); exit; end;
 for i:=3 to n do
  begin
   c:=last[i];
   square_sums[i]:=square_sums[c]+i/sum*i;
   if c<>0 then
    square_sumt:=square_sumt+i/sum*i*f[c]+square_sums[c]-2*i/sum*sums[c];
  end;
 ave:=sumt/sum;
 ans:=ave*ave-2*ave*sumt/sum+square_sumt;
 writeln(ans:0:6);
end;

begin
 assign(input,'tri.in');
 assign(output,'tri.out');
 reset(input);
 rewrite(output);
 readln(t);
 while t>0 do
  begin
   main;
   dec(t);
  end;
 close(input);
 close(output);
end.