AC

pe vol31 151~155

ben1222 — Sun, 01 May 2011 13:34:54 +0800

最近好像没什么干劲啊…

Problem 151 : C++,裁纸,数学期望
Problem 152 : C++,1/2用整数的平方的倒数的和来表示的方法数
Problem 153 : C++,整除n的高斯整数之和
Problem 154 : C++,(x+y+z)²⁰⁰⁰⁰⁰ 中有多少个系数是10¹²的整数倍
Problem 155 : C++,用18个相同容量的电容串联并联等能组成多少种不同的电容值

下面是详细内容:

Problem 151
Paper sheets of standard sizes: an expected-value problem.
Project Euler 151
一开始只有一张A1的纸,每次随机取出一张纸
如果不是A5则按题意切,切完拿走一张A5,其余放回
如果是A5则拿走
除去第一次和最后一次,拿的时候发现只剩一张纸的次数的数学期望是多少

每个状态算一遍就行了
对每一天来说,记录 {各大小纸张的数量,至今为止遇到只剩一张纸的情况的次数} = 出现概率
状态不多,用map就够用了
从第一天开始依次计算下一天的所有状态的出现概率
然后(概率*次数)求和即可

Problem 152
Writing 1/2 as a sum of inverse squares
把1/2表示成许多不同整数的平方的倒数的和
求从2~80这几个整数取数,按上述计算,有多少种方法表示1/2

既然最后要表示成1/2,那么分母中除了2以外的质数因子都要能消除才行
怎么才能消除呢…
1/a+1/b=(a+b)/(a*b)
如果a和b互质的话,(a+b)和(a*b)也互质
c/a+d/b=(b*c+a*d)/(a*b)
如果a和b互质且c/a和d/b是既约分数的话,(b*c+a*d)和(a*b)也互质…

因此为了消除分母上的某个质因子,这个质因子至少要在两项中出现
分母只能取2~80的平方,超过40的质数最多只有一个…所以这些数就不考虑了
而超过7的质数之间又是不相关的(2~80中没有哪个数是两个超过7的质数的积)
既然是不相关的就可以单独拿出来看看…能有哪几种消除该质数因子的方法

写段程序检测下来发现含有大于7的质数的因子的数中只有这么一组可以消除该质数的:

	1/144=1/169+1/1521+1/2704=1/13²+1/39²+1/52²

那么还剩下39个数:

2,3,4,5,6,7,8,9,10,12,14,15,16,18,20,21,24,25,27,28,
30,32,35,36,40,42,45,48,49,50,54,56,60,63,64,70,72,75,80

算上1/144这个,还剩40个数,其中12出现2次

既然这结果这么好,那么带7这个因子的数也拿来组合一下试试…
这下就有好多种消除质因子7的组合了,有19个
剩下的数字有30个:

2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,
30,32,36,40,45,48,50,54,60,64,72,75,80,12

本来以为这19个组合里面如果两两没有交集的话,还要考虑选取这2个的情况
实际上多虑了,因为两个都选取的情况已经包含在这些组合里面了

那么在这19个数任取一个或者一个都不取…然后和带有5的因子的数组合起来看看消除5这个因子有多少种组合
有475个组合,其中不相同的数有252个
还剩18个数字:
2,3,4,6,8,9,12,16,18,24,27,32,36,48,54,64,72,12

对于3和2的因子也如法炮制,对有重复的数字做个乘法就行了,不要重新算
结果好小…答案正确

Problem 153
Investigating Gaussian Integers
定义高斯整数为实部和虚部都是整数的复数
定义若高斯整数a和高斯整数b满足b/a也是高斯整数,那么称a能整除b
s(n)表示所有能整除n的实部为正的高斯整数的和(即满足n/x是高斯整数的x的和)
求sum{s(n) ,n=1..10^8}

–下面的除法是整数除法
设s(n)=sr(n)+si(n),其中sr(n)表示整数部分的和,si(n)表示虚部不为0的复数部分的和
整数部分的话,很简单

sum{sr(n),n=1..N}=sum{(N/n)*n,n=1..N}

可以优化一下用O(sqrt(N))的复杂度算出来

不过复数不太好办
设最小的能被a+bi整除的正整数是m(a,b)

m(a,b)/(a+bi)=m(a,b)*(a-bi)/(a²+b²)
因此
m(a,b)=(a²+b²)/gcd(a,b)

于是能被a+bi整除的正整数必须是m(a,b)的倍数
共轭复数是成对出现的,即a+bi可以整除n的话a-bi也可以,而加起来正好是2*a,虚部被抵消了
因此只计算b>0的情况就行了

si(n)=2*sum{a| n%m(a,b)==0, b>0}

分成a=b,a<b和a>b的的情况讨论得到

si(n)=2*sum{a| n%(2*a)==0}+2*sum{a+b| n%m(a,b)==0,a>b>0}

求和后就成了

sum{si(n),n=1..N}=2*sum{(N/(2*a))*a| a=1..N}+2*sum{(N/m(a,b))*(a+b)| 0<b<a≤N}

前半部分和之前算整数部分之和的方法差不多,后半部分怎么办呢
可以开个10^8大的数组,r[i]记录N/i的系数
对每个0<b<a<sqrt(N) ,gcd(a,b)=1的(a,b),给r[(a*a+b*b)*k]+=(a+b)*2*k
相信只有O(N)个这样的项需要处理,gcd复杂度算是O(logN)那么总复杂度是O(N*logN)

感觉内存用的太多,空间复杂度应该可以简化到O(sqrt(N))级别的
r0[1..sqrt(N)]记录N/i的系数
r1[1..sqrt(N)]记录N/k=i的系数

内存消耗不到1MB,耗时20秒
答案正确

Problem 154
Exploring Pascal’s pyramid.
金字塔形堆放的球,顶端的球数字为1,其他球的数字是它上面相邻的3个球的数字之和
也就是每个球上的数字是三项式(x+y+z)ⁿ展开后的系数
问(x+y+z)²⁰⁰⁰⁰⁰ 中有多少个系数是10¹²的整数倍

让人想起Problem 148…不过又是3维又不是质数…
先展开多项式吧
设w=x+y,那么

(x+y+z)ⁿ=(w+z)ⁿ=sum{C(n,i)*wⁱ*z^(n-i)}
	=sum{C(n,i)*(x+y)ⁱ*z^(n-i)}
	=sum{C(n,i)*C(i,j)*x^j*y^(i-j)*z^(n-i)}

设a=j,b=i-j,c=n-i,分别表示x,y,z的幂次

C(n,i)*C(i,j)=n!*i!/(i!*(n-i)!*j!*(i-j)!)
	=n!/(a!*b!*c!)

要求这个系数是10¹²的整数倍啊…
就是看这个数中2的幂次和5的幂次哪个低咯…低的那个达到12就ok
1~n每个数的阶乘的2的幂次和5的幂次都可以预处理下
枚举a和b,c=n-a-b,复杂度是O(n^2)…虽然不是不能接受…但是不是有点慢…

尝试从小数据中找规律…结果失败了…

那就看看能不能做些常数优化…
设a<=b<=c,减少计算次数
枚举范围是a=[0,n/3],b=[a,(n-a)/2]
那么a,b,c都不相等的时候,n!/(a!*b!*c!)值共出现了6次
a=b或者b=c的时候,出现3次
a=b=c的时候,出现1次

34秒,答案正确

Problem 155
Counting Capacitor Circuits.
设D(n)表示最多用n个相同容量的电容串联并联等能组成多少种不同的电容值
求D(18)

写个处理分数的类…然后暴力
设c(n)表示用n个相同容量的电容可以组成哪些电容值
那么枚举c(i),c(n-1)两两配对进行串联或并联,然后排序去重…
48秒
答案正确

大数取余(二) 剩余定理

ben1222 — Thu, 07 Apr 2011 21:06:50 +0800

上次的文章中主要说的都是各种运算对应的取余的计算方式
但是没有提到关于模数的优化

计算a%m时,如果m能小一些的话,通常计算起来就会方便很多
比如可以用更小的内存进行预处理,或者找到循环节之类的很好的东西

在做了Project Euler 330后深刻的感受到了这点…
题中要求对77777777取模,可以看到这和平时常用的对某个质数取模不同,这回是个合数
77777777 = 7 * 11 * 73 * 101 * 137
在过去,看到要对合数取模感觉就不太爽,因为有些性质只有质数才能用,比如组合数取模

这题搞了很久没有进展才想起来还有个剩余定理
对于两两互质的一组整数{m₁,m₂,…,m_n}
设m=m₁*m₂*…*m_n, 也就是m是这组数的乘积
设b_i=a%m_i, 也就是b_i是a对m_i取余的余数
已知{b₁,b₂,…,b_n}的话,可以求得a%m
具体做法是找到c_i=(m/m_i)*k使得c_i%m_i==1
因为两两互质,所以k只要枚举1~m_i-1就一定能找到满足条件的c_i
最后a%m=sum{c_i*b_i}%m

因此当m可以分解质因数时可以分成两两互质的一组数,然后通过计算b_i来得到a%m
这样做有2个好处
一个是有可能得到m_i是质数,可以应用质数所特有的方法
另一个是模数小了,这样如果有循环节的话就更容易发现,这就更爽了

define型计算几何模板

ben1222 — Wed, 23 Feb 2011 15:25:48 +0800

本想在寒假期间写点什么的…无奈域名出了点问题…
导致这一个半月里面网站一直处于无法访问的状态…

想起以前积累的计算几何模板…想差不多可以整理下发出来…
不过充斥着define…以及define的奇怪用法…可以说完全是C语言风格的代码…
总之这是一个基于数组和define而不是类和template的模板…会有人看吗…

如果每段代码都注解出一步步怎么得到的话…工程浩大…
结果心血来潮…给几乎每段代码画了个图片…
蓝色是已知,红色是所求…
画了几个才发现…有些还真难用图片表达啊…

那么就开始吧…

先是关于点,线,圆的形式的说明

说明:
point,p,vector,v,点,向量:
  array[2]  [0]:x,dx  [1]:y,dy
line,vl,直线:
  array[3]  [0]:a  [1]:b  [2]:c
  对应直线方程ax+by+c=0
circle,o,圆:
  array[3]  [0]:x  [1]:y  [2]:r

常量:圆周率和微小值

#define PI 3.14159265358979
#define EPS 1e-10

然后是向量运算

#define IM(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1])
#define IM2(ps1,pe1,ps2,pe2) ((pe1[0]-ps1[0])*(pe2[0]-ps2[0])+(pe1[1]-ps1[1])*(pe2[1]-ps2[1]))
/* 外积大于0,则v2在v1左手方向 */
#define OM(v1,v2) (v1[0]*v2[1]-v1[1]*v2[0])
#define OM2(ps1,pe1,ps2,pe2) ((pe1[0]-ps1[0])*(pe2[1]-ps2[1])-(pe1[1]-ps1[1])*(pe2[0]-ps2[0]))
#define SU(v1,v2,v3) v3[0]=v2[0]-v1[0];v3[1]=v2[1]-v1[1];

距离:

/*两点间距离平方*/
#define DIS2(p1,p2,ret) {double __v[2];SU(p1,p2,__v) ret=IM(__v,__v);}
/*两点间距离*/
double dis(long*p1,long*p2){
	double v[2];
	SU(p1,p2,v);
	return sqrt(IM(v,v));
}

关于圆:
三点求圆心:


/* 三点求圆心 */
#define SQS(n,d) n=(d)*(d)
#define ARCC(x1,y1,x2,y2,x3,y3,x0,y0) {\
	double SQS(__y1,y1),SQS(__y2,y2),SQS(__y3,y3),SQS(__x1,x1),SQS(__x2,x2),SQS(__x3,x3); \
	x0=((y3-y1)*(__y2-__y1+__x2-__x1)+(y2-y1)*(__y1-__y3+__x1-__x3))*0.5/ \
		(double)((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)); \
	y0=((x3-x1)*(__x2-__x1+__y2-__y1)+(x2-x1)*(__x1-__x3+__y1-__y3))*0.5/ \
		(double)((y2-y1)*(x3-x1)-(y3-y1)*(x2-x1));}
#define ARCCP(p1,p2,p3,pc) ARCC(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],pc[0],pc[1])

两点和直径求圆心:


/*
给定圆上两点和圆的直径求2个候选的圆心位置,不判断两点间距离是否超过直径
使用:math.h
*/
#define PDCC(x1,y1,x2,y2,d,rx1,ry1,rx2,ry2) {\
	double SQS(_x,x2-x1),SQS(_y,y2-y1);\
	double _u=sqrt(d*d/(_x+_y)-1);\
	rx1=((x1+x2)+(y2-y1)*_u)/2;\
	ry1=((y1+y2)-(x2-x1)*_u)/2;\
	rx2=((x1+x2)-(y2-y1)*_u)/2;\
	ry2=((y1+y2)+(x2-x1)*_u)/2;}
#define PDCCP(p1,p2,d,r1,r2) PDCC(p1[0],p1[1],p2[0],p2[1],d,r1[0],r1[1],r2[0],r2[1])

圆外一点到圆的切点:


/*
求圆外一点到圆的2条切线的2个切点
po是圆心坐标,pp是圆外的一点,r是圆的半径,ret1和ret2存放2个切点坐标
不检测点是否在圆内
使用:向量运算,math.h
*/
#define TPC(p1,op,p2) (po[p1]+r*(r*d1[p1] op d1[p2]*sd)/s1)
void tanpoint(double*po,double*pp,double r,double*ret1,double*ret2){
	double d1[2],s1,sd;
	SU(po,pp,d1)
	s1=IM(d1,d1);
	sd=s1-r*r;
	if(sd<0)sd=-sd;
	sd=sqrt(sd);
	ret1[0]=TPC(0,+,1);
	ret1[1]=TPC(1,-,0);
	ret2[0]=TPC(0,-,1);
	ret2[1]=TPC(1,+,0);
}

圆相交面积:


/*
圆相交面积
使用:向量运算,距离,常量,math.h
*/
double cosine(double b1, double b2, double b3){
	double b = b2 * b2 + b3 * b3 - b1 * b1;
	return b / b2 / b3 / 2;
}
double circle_ins(double*o1,double*o2){
	double dc=dis(o1,o2),r2=(o1[2]>o2[2]?o1[2]:o2[2]),r1=(o1[2]<o2[2]?o1[2]:o2[2]);
	if(r2-r1>=dc)return PI*r1*r1;
	else if(dc>=r2+r1)return 0;
	else{
		double a1,a2,ang1,ang2,ret=0,s1=PI*r1*r1,s2=PI*r2*r2;
		a1=cosine(r2,r1,dc);
		a2=cosine(r1,r2,dc);
		ang1=acos(a1);
		ang2=acos(a2);
		if(a1>0)ret=(r1*r1*ang1-r1*r1*sin(ang1*2)/2);
		else ret=(r1*r1*ang1+r1*r1*sin((PI-ang1)*2)/2);
		if(a2>0)ret+=(r2*r2*ang2-r2*r2*sin(ang2*2)/2);
		else ret+=(r2*r2*ang2+r2*r2*sin((PI-ang2)*2)/2);
		return ret;
	}
}

关于直线:
直线解析式:


/*
给定两点,求其所在直线的解析式
ret[0]=a,ret[1]=b,ret[2]=c
ax+by+c=0
使用:向量运算
*/
#define PT_LINE(p1,p2,ret) {ret[1]=p2[0]-p1[0];ret[0]=p1[1]-p2[1];ret[2]=-IM(p1,ret);}

点关于直线的对称点:


/*
点关于直线的对称点
vl:直线解析式,pt:点,ret:返回对称点
使用:向量运算
*/
#define PL_SYM(vl,pt,ret) {\
	double __d=IM(vl,vl); \
	ret[0]=((vl[1]*vl[1]-vl[0]*vl[0])*pt[0]-2*vl[0]*(vl[1]*pt[1]+vl[2]))/__d; \
	ret[1]=((vl[0]*vl[0]-vl[1]*vl[1])*pt[1]-2*vl[1]*(vl[0]*pt[0]+vl[2]))/__d;}

直线交点:


/*
判断两条直线[ps1,pe1],[ps2,pe2]的交点情况:
	0:平行 1:相交 2:重叠
如果有交点,返回值在ret中
使用:向量运算
*/
long lineins(long*ps1,long*pe1,long*ps2,long*pe2,double*ret){
	long v1[2],v2[2];
	v1[0]=pe1[0]-ps1[0];v1[1]=pe1[1]-ps1[1];
	v2[0]=pe2[0]-ps2[0];v2[1]=pe2[1]-ps2[1];
	long mul=OM(v1,v2);
	if(mul==0){
		v2[0]=pe2[0]-ps1[0];
		v2[1]=pe2[1]-ps1[1];
		return (OM(v1,v2)==0)?2:0;
	}else{
		double t=((ps2[0]-ps1[0])*v2[1]+(ps1[1]-ps2[1])*v2[0])/(double)mul;
		ret[0]=ps1[0]+v1[0]*t;
		ret[1]=ps1[1]+v1[1]*t;
		return 1;
	}
}

点到直线的距离:


/*
点pt到直线[pl1,pl2]的距离平方
使用:向量运算
*/
double dis2pl(double*pt,double*pl1,double*pl2){
	double v1[2],v2[2];
	SU(pt,pl1,v1)SU(pl1,pl2,v2)
	double om=OM(v1,v2);
	return om*om/IM(v2,v2);
}

直线和圆的交点:


/*
直线[pl1,pl2]与圆[po,r]的交点,返回值为交点个数
ax+by+c=0;
(x-ox)^2+(y-oy)^2=r^2
x=( b(b*ox-a*oy)-a*c土b*sqrt(r^2-dis2pl()^2) )/( a^2+b^2 )
y=( a(a*oy-b*ox)-b*c干a*sqrt(r^2-dis2pl()^2) )/( a^2+b^2 )
使用:向量运算,常量,直线解析式,math.h
*/
long lcins(double*po,double r,double*pl1,double*pl2,double*ret1,double*ret2){
	double vl[3],tmp1,tmp2;
	PT_LINE(pl1,pl2,vl)
	tmp1=IM(vl,po)+vl[2];
	double ab=IM(vl,vl),mdis2=ab*r*r-tmp1*tmp1;
	if(mdis2<=-EPS)return 0;
	else{
		tmp1=(vl[1]*OM(po,vl)-vl[0]*vl[2]);
		tmp2=(vl[0]*OM(vl,po)-vl[1]*vl[2]);
		if(mdis2<EPS){
			ret1[0]=tmp1/ab;ret1[1]=tmp2/ab;
			return 1;
		}else{
			double sq=sqrt(mdis2);
			ret1[0]=(tmp1+vl[1]*sq)/ab;ret1[1]=(tmp2-vl[0]*sq)/ab;
			ret2[0]=(tmp1-vl[1]*sq)/ab;ret2[1]=(tmp2+vl[0]*sq)/ab;
			return 2;
		}
	}
}

直线和椭圆的交点:


/*
直线[pl1,pl2]与椭圆的交点,返回值为交点个数
ax+by+c=0
(x/d)^2+(y/e)^2=1
x=( -c*ad*d土be*d*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
y=( -c*be*e干ad*e*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
使用:向量运算,常量,直线解析式,math.h
*/
long leins(double d,double e,double*pl1,double*pl2,double*ret1,double*ret2){
	double vl[3];
	PT_LINE(pl1,pl2,vl)
	double ab=vl[0]*vl[0]*d*d+vl[1]*vl[1]*e*e,delta=ab-vl[2]*vl[2];
	if(delta<=-EPS)return 0;
	else if(delta<EPS){
		ret1[0]=-vl[2]*vl[0]*d*d/ab;
		ret1[1]=-vl[2]*vl[1]*e*e/ab;
		return 1;
	}else{
		double sq=sqrt(delta);
		ret1[0]=(-vl[2]*vl[0]*d*d+vl[1]*e*d*sq)/ab;
		ret1[1]=(-vl[2]*vl[1]*e*e-vl[0]*e*d*sq)/ab;
		ret2[0]=(-vl[2]*vl[0]*d*d-vl[1]*e*d*sq)/ab;
		ret2[1]=(-vl[2]*vl[1]*e*e+vl[0]*e*d*sq)/ab;
		return 2;
	}
}

关于线段:

/*
整数点间连一条线,问穿过多少个单位正方形:
	=a+b-gcd(a,b)
扩展到三维:
	=a+b+c-gcd(a,b)-gcd(a,c)-gcd(b,c)+gcd(a,b,c)
*/

判断线段是否相交:


/*
判断两线段是否相交,跨立实验中用>则不包括端点,用>=则包括端点
使用:向量运算
*/
#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MN(a,b,t) (MAX(ps##a[t],pe##a[t])>=MIN(ps##b[t],pe##b[t]))
long segins(long*ps1,long*pe1,long*ps2,long*pe2){
	return MN(1,2,0)&&MN(2,1,0)&&MN(1,2,1)&&MN(2,1,1) && //快速排斥实验
		(OM2(ps1,ps2,ps1,pe1)*OM2(ps1,pe1,ps1,pe2)>=0) &&
		(OM2(ps2,ps1,ps2,pe2)*OM2(ps2,pe2,ps2,pe1)>=0); //跨立实验
}

点到线段的垂足:


/*
点pt到线段[ps,pe]所在直线的垂足,返回true表示垂足在线段上,垂足坐标存在ret中
使用:向量运算
*/
long plfoot(double*pt,double*ps,double*pe,double*ret){
	double vt[2],t,ret[2];
	SU(ps,pe,vt)
	t=(-(ps[0]-pt[0])*(pe[0]-ps[0])-(ps[1]-pt[1])*(pe[1]-ps[1]))/IM(vt,vt);
	ret[0]=ps[0]+t*(pe[0]-ps[0]);
	ret[1]=ps[1]+t*(pe[1]-ps[1]);
	return !((ret[0]<ps[0]&&ret[0]<pe[0])||(ret[0]>ps[0]&&ret[0]>pe[0])||
		(ret[1]<ps[1]&&ret[1]<pe[1])||(ret[1]>ps[1]&&ret[1]>pe[1]));
}

点到线段的距离:


/*
点pt到线段[pl1,pl2]的距离的平方
使用:向量运算,常量,距离,点到直线距离
*/
double dis2ps(double*pt,double*pl1,double*pl2){
	double ret;
	if(IM2(pl1,pt,pl1,pl2)<-EPS)DIS2(pt,pl1,ret)
	else if(IM2(pl2,pt,pl2,pl1)<-EPS)DIS2(pt,pl2,ret)
	else ret=dis2pl(pt,pl1,pl2);
	return ret;
}

角平分线上一点:


/*
求2个线段[pc,pl]和[pc,pr]组成的夹角的角平分线上的一个点
使用:向量运算,math.h
*/
void hangle(double*pl,double*pc,double*pr,double*ret){
	double vl[2],vr[2];
	SU(pc,pl,vl)
	SU(pc,pr,vr)
	double dl=sqrt(IM(vl,vl)),dr=sqrt(IM(vr,vr));
	vl[0]/=dl;vl[1]/=dl;
	vr[0]/=dr;vr[1]/=dr;
	ret[0]=(vl[0]+vr[0])/2.0+pc[0];
	ret[1]=(vl[1]+vr[1])/2.0+pc[1];
}

关于三角形:
三角形面积:


/*
三角形面积之海伦公式:
p=(a+b+c)/2
s=sqrt((p-a)*(p-b)*(p-c)*p)
*/

三角形内切圆半径:


/*
三角形内切圆半径:
直角三角形:r=(a+b-c)/2,其中a,b是直角边长,c是斜边长
一般三角形:r=2S/(a+b+c),其中S是三角形面积
*/

三角形内周长一定时的最大面积:


/*
三角形内周长一定的最大面积
triple=x+y+z;
half=triple*0.5;
area=sqrt(half*(half-x)*(half-y)*(half-z));
R=area*2.0/triple;
if(x+y+z<=role)
	max=area;
else if(2.0*PI*R>=role)
	max=role*role/(4.0*PI);
else
{
	r=(triple-role)/(triple/R-2.0*PI);
	max=area+PI*r*r-(r*r*area/(R*R));
}

*/

三角形重心:

/* 三角形重心 */
#define CGT(p1,p2,p3,pc) {pc[0]=(p1[0]+p2[0]+p3[0])/3;pc[1]=(p1[1]+p2[1]+p3[1])/3;}

判断点是否在三角形内:


/*
外积判断点是否在三角形内
在三角形内的返回2,边上的1,外面的0
使用:向量运算
*/
long pint(long*p1,long*p2,long*p3,long*pm){
	long v1[2],v2[2],t,l;
	SU(p1,p2,v1)SU(p1,pm,v2)
	if((t=OM(v1,v2))==0)return 1;
	t=t>0?1:-1;
	SU(p2,p3,v1)SU(p2,pm,v2)
	l=t*OM(v1,v2);
	if(l==0)return 1;else if(l<0)return 0;
	SU(p3,p1,v1)SU(p3,pm,v2)
	l=t*OM(v1,v2);
	if(l==0)return 1;else if(l<0)return 0;
	return 2;
}

三角形内一点到3边的距离的平方和最小:(这个图片我实在想不出怎么画…)

/*
三角形内一点到3边的距离的平方和最小
void (double pts[3][2])
	double x,y,vl[3][3],d2[3];
	long i;
	for(i=0;i<3;i++){
		PT_LINE(pts[i],pts[(i+1)%3],vl[i]);
		d2[i]=IM(vl[i],vl[i]);
	}
	double X=0,Y=0,C=0,U=0,B=0;
	for(i=0;i<3;i++){
		X+=vl[i][0]*vl[i][1]/d2[i];
		Y+=vl[i][1]*vl[i][1]/d2[i];
		C+=vl[i][1]*vl[i][2]/d2[i];
	}
	for(i=0;i<3;i++){
		double t=vl[i][0]*Y-vl[i][1]*X;
		U+=t*(vl[i][2]*Y-vl[i][1]*C)/d2[i];
		B+=t*t/d2[i];
	}
	x=-U/B;
	y=-(X*x+C)/Y;
*/

判断矩形能否斜放在另一个矩形内:


/*
判断矩形能否斜放在另一个矩形内
minh=(2xya+d*(x*x-y*y))/(x*x+y*y)  (x>y)
d=sqrt(x*x+y*y-a*a)
return minh<=b

不判断面积大小和快速排斥
使用:math.h
*/
#define I __int64
long rcinrc(long w1,long h1,long w2,long h2){
	if((w1<h1?w1:h1)<(w2<h2?w2:h2))return 0;
	I sq=(I)w2*(I)w2+(I)h2*(I)h2,l=2*(I)w2*(I)h2*(I)w1,a2=(I)w1*(I)w1;
	if(w2==h2||sq==a2)return l<=sq*h1;
	else{
		if(w2<h2){long t=w2;w2=h2;h2=t;}
		if(sq<a2)return 0;
		I il=((I)w2*(I)w2-(I)h2*(I)h2);
		double dl=l+sqrt((double)(sq-a2))*il,dr=sq*h1;
		return dl<=dr;
	}
}

三角形和圆的公共部分面积:


/*
点pt在线段[pl1,pl2]的位置
IM2(pl1,pt,pl1,pl2)<-EPS  点在线段外,pl1一侧
IM2(pl2,pt,pl2,pl1)<-EPS  点在线段外,pl2一侧
PTINSEG:点在线段所在直线的前提上,判断点是否在线段内,线段加长一段微小量
PTINSEGS:同上,线段缩小一段微小量
使用:向量运算,常量
*/
#define PTINSEG(pt,pl1,pl2) ((IM2(pl1,pt,pl1,pl2)>=-EPS)&&(IM2(pl2,pt,pl2,pl1)>=-EPS))
#define PTINSEGS(pt,pl1,pl2) ((IM2(pl1,pt,pl1,pl2)>EPS)&&(IM2(pl2,pt,pl2,pl1)>EPS))

/*
2点和圆心组成的三角形与这个圆的公共部分面积
r是圆半径,原点是圆心
ps和pe是三角形另外两个点,逆时针为正
使用:向量运算,常量,直线和圆的交点,math.h
*/
#define VRAD(v,ret) {if(v[0]==0){if(v[1]>0)ret=PI/2;else ret=PI*3/2;}\
	else{double _t=atan(v[1]/v[0]);if(v[0]<0)_t=PI+_t;ret=(_t<0?_t+PI*2:_t);}}
#define MARC(rad) {if(rad>PI)rad-=PI*2;if(rad<-PI)rad+=PI*2;}
double area_c2p(double r,double*ps,double*pe){
	double d1=IM(ps,ps),d2=IM(pe,pe),r2=r*r;
	if(d1<=r2&&d2<=r2)return OM(ps,pe)/2;
	double rad1,rad2,po[2],ins1[2],ins2[2],drad;
	VRAD(ps,rad1)VRAD(pe,rad2)
	drad=rad2-rad1;MARC(drad)
	po[0]=po[1]=0;
	long nins=lcins(po,r,ps,pe,ins1,ins2);
	if(nins<=1)return drad/2*r*r;
	if((IM2(ps,ins1,ps,pe)<-EPS&&IM2(ps,ins2,ps,pe)<-EPS)||\
		(IM2(pe,ins1,pe,ps)<-EPS&&IM2(pe,ins2,pe,ps)<-EPS))return drad/2*r*r;
	else if(PTINSEG(ins1,ps,pe)&&PTINSEG(ins2,ps,pe)){
		double radi1,radi2,drad2,ret;
		if(PTINSEGS(ins1,ps,ins2)){VRAD(ins1,radi1)VRAD(ins2,radi2)ret=OM(ins1,ins2)/2;}
		else{VRAD(ins2,radi1)VRAD(ins1,radi2)ret=OM(ins2,ins1)/2;}
		drad=radi1-rad1;MARC(drad)
		drad2=rad2-radi2;MARC(drad2)
		return ret+(drad+drad2)/2*r*r;
	}else if(PTINSEG(ps,ins1,ins2)&&PTINSEG(pe,ins1,ins2))return OM(ps,pe)/2;
	else{
		double*pins,*pout,radi;
		if(IM2(ps,ins1,ps,pe)<-EPS||IM2(pe,ins1,pe,ps)<-EPS){pins=ins2;pout=ins1;}
		else{pins=ins1;pout=ins2;}
		VRAD(pins,radi)
		if(IM2(pe,pout,pe,ps)<-EPS){
			drad=radi-rad1;MARC(drad)
			return OM(pins,pe)/2+drad/2*r*r;
		}else{
			drad=rad2-radi;MARC(drad)
			return OM(ps,pins)/2+drad/2*r*r;
		}
	}
}

凸包:


/*排序用*/
#define POINT struct TPoint
POINT{long xy[2];}pstart;
long operator<(POINT a,POINT b){
	long t=OM2(pstart.xy,a.xy,pstart.xy,b.xy);
	if(t==0){
		long d1,d2;
		DIS2(pstart.xy,a.xy,d1)
		DIS2(pstart.xy,b.xy,d2)
		return d1<d2;
	}else return t>0;
}
/*
O(nlogn)的凸包,返回值是凸多边形的顶点数,存在pts中,逆时针方向
使用:向量运算,距离,std::sort(algorithm)
*/
long convex(POINT*pts,long num,long*stack){
	long i,minp,minx=0x7fffffff,miny;
	for(i=0;i<num;i++){
		if(pts[i].xy[0]<minx){minx=pts[i].xy[0];miny=pts[i].xy[1];minp=i;}
		else if(pts[i].xy[0]==minx&&pts[i].xy[1]<miny)miny=pts[minp=i].xy[1];
	}
	if(minx==0x7fffffff)return 0;
	POINT _t=pts[minp];
	pts[minp]=pts[0];
	pts[0]=_t;
	pstart=pts[0];
	sort(pts+1,pts+num);
	long top=3;
	for(i=0;i<3;i++)stack[i+1]=i;
	for(i=3;i<num;i++){
		while(OM2(pts[stack[top-1]].xy,pts[stack[top]].xy,pts[stack[top-1]].xy,pts[i].xy)<0)top--;
		stack[++top]=i;
	}
	for(i=1;i<=top;i++)pts[i-1]=pts[stack[i]];
	return top;
}

凸洞:见这里

pe vol30 146~150

ben1222 — Fri, 29 Oct 2010 10:16:44 +0800

最近好迷茫…我想做什么?我能做什么?
没有排名…没有考试…没有比较…不知道自己几斤几两…该往什么方向努力…
有能测出人各方面的长处短处的东西会怎么样呢?这样设想着…

Project Euler
pe第三页拖拖拉拉的终于切完了
(最近的新题貌似不太给力啊…老题还是有不少有意思的题目…)

Problem 146 : C++,n^2+k连续素数,减少枚举量
Problem 147 : C++,数矩形,推出式子后就是简单的dp了
Problem 148 : C++,杨辉三角与分形,图形很美,画出图形后算起来就简单多了
Problem 149 : C++,矩阵中行、列、对角线方向的最大子段和
Problem 150 : C++,最大子三角形和,各行预处理下即可

下面是详细内容:

Problem 146
Investigating a Prime Pattern
Project Euler 146
找n使其满足n^2+1,n^2+3,n^2+7,n^2+9,n^2+13,n^2+27是连续的质数
最小的正整数n是10
求n<150,000,000中满足上面条件的n的和

首先,n肯定是偶数
而且不能有3,7,13 的因子…
然后…枚举判断试试看
速度好慢…不过这样的数好稀少…一百万中才3个…
10
315410
927070

貌似必须是10的倍数?
如果末尾是2的话,n^2+1的末尾就是5了,不是素数
如果末尾是4的话,n^2+9的末尾就是5了
末尾是6的话,n^2+9的末尾也是5
末尾是8的话,n^2+1的末尾也是5
所以n必须是10的倍数
因此n^2+5,n^2+15,n^2+25都一定是合数了,不用判断他们是否是合数
速度加快了近4倍…

n^2末尾2位是0了…
根据3的倍数的各位和是3的倍数这一点可以知道…n^2各位和必须是3k+1,否则加上末尾的1,3,7,9,13,27的各位和就会是3的倍数了
因此,n^2+11,n^2+17,n^2+23必定能被3整除,不用判断了

跑了1.5小时才出来结果…答案正确
貌似还可以用n对7的余数来进一步减少判断量

Problem 147
Rectangles in cross-hatched grids
a*b的网格,每个单位方格中有对角线
这样的网格中有多少个矩形呢?包括斜过来的
设这个函数是f(a,b)
求sum{f(a,b),a<=47,b<=43},(a,b是正整数)

想DP…可是不好办…因为那个斜的…
首先,f(a,b)=f(b,a)这是肯定的
那么可以想象a*b的网格由a*(b-1)的网格2个,中间重叠了a*(b-2)的部分
f(a,b)=f(a,b-1)*2-f(a,b-2)+g(a,b)
这里的g(a,b)是那些同时占据到2条a*1的那些矩形,这样的矩形数量好算些…吧
增加的平放的矩形很简单,是a*(a+1)/2个
增加的斜放的矩形嘛…设为h(a,b)个
由于每条只有1的宽度,所以应该也不难算
看看有哪些是有可能组成矩形的

下侧的对称

中间2种情况其实是对称的,算作一种
从左到右依次编号(1),(2),(2),(3)
那么看看上侧和下侧的组合的情况:h(a,b)=h_11(a,b)+h_12(a,b)+h_13(a,b)+h_21(a,b)+..+h_33(a,b)

(1)-(1):

    h_11(n,2)=n
    h_11(n,3)=(n-1)*2+(n-2)   (n>=2)
    h_11(n,4)=(n-2)*2+(n-3)*2+(n-2)  (n>=2)
    h_11(n,5)=(n-3)*2+(n-4)*2+(n-3)*2+(n-4)  (n>=4)
    h_11(a,b)=
              {(a*2-(b*2-3))*(b-2)+(a-(b-2))  b is even
              {(a*2-(b*2-3))*(b-1)-(a-(b-1))  b is odd
             =(a*2-b*2+3)*(b-2)+(a-b+2)

(1)-(2), (2)-(1):限制住必须是宽为sqrt(2)/2的矩形了

    h_12(a,b)=h_21(a,b)=h_11_0(a-1,b)=((a-1)-(b-2))*2
        =(a-b+1)*2

(1)-(3), (3)-(1):必须宽超过sqrt(2)/2的矩形,和上面(1)-(2),(2)-(1)那种合并起来看就比较方便

    h_12(a,b)+h_13(a,b)=h_21(a,b)+h_31(a,b)
        =(a-(b-1))*2*(b-1)
        =(a-b+1)*2*(b-1)

(2)-(3), (3)-(2):这个不存在满足的情况
(2)-(2), (3)-(3):这2个也合在一起看比较方便

    h_22(a,b)+h_33(a,b)
    b=4:(a-b)*2+(a-b+1)*2+(a-b)*2+(a-b+1)
    b=5:(a-b)*2+(a-b+1)*2+(a-b)*2+(a-b+1)*2+(a-b)
    h_22(a,b)+h_33(a,b)=(a-b)*b+(a-b+1)*(b-1)

连在一起就是:

h(a,b)=h_11(a,b)+(h_12(a,b)+h_13(a,b))*2+(h_22(a,b)+h_33(a,b))
    =(a*2-b*2+3)*(b-2)+(a-b+2)+(a-b+1)*2*(b-1)*2+(a-b)*b+(a-b+1)*(b-1)
    =8*a*b-8*a-8*b^2+16*b-9
g(a,b)=a*(a+1)/2+8*a*b-8*a-8*b^2+16*b-9

然后就是简单的dp了,答案正确

Problem 148
Exploring Pascal’s triangle.
杨辉三角中前10^9行中不能被7整除的数有多少

杨辉三角中所有数都用它除7的余数来表示的话,就是找不为0的项的个数
设A={
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 3 3 5 1
1 6 1 6 1 6 1
}
也就是模7的杨辉三角前7行,有28个非7倍数的数
那么前14行是
A
A A
中间空余部分倒三角形的区域全是0,也就是是7的倍数
前49行是
B={
A
A A
A 2A A
A 3A 3A A
A 4A 6A 4A A
A 5A 3A 3A 5A A
A 6A 1A 6A 1A 6A A
}
由于7是质数,所以n*A(n<7)不会出现7的倍数
来张图片好理解一些:很像Sierpinski三角形吧

其中灰色的部分是7的倍数

算非7倍数的项很方便:
前7^2行有28^2个非7倍数的项
前7^3行有28^3个非7倍数的项
以此类推

剩下的就简单了,O(logn)的复杂度搞定吧
瞬出,答案正确

Problem 149
Searching for a maximum-sum subsequence.
2000*2000的矩阵(有负数)中找出行、列、对角线方向的连续项之和的最大值

对每整行、整列、整对角线分别做最大子段和
复杂度O(n^2)
不到一秒就出结果,答案正确

Problem 150
Searching a triangular array for a sub-triangle having minimum-sum.
堆成一千行三角形的数字,找一个子三角形,使得其包含的数字之和最小

枚举子三角形的上顶点,然后往下扩展
每行都预处理下,以便快速得到该行的各区间和
O(n^3),隐藏的常系数较小
5秒搞定,答案正确

flash Electric box 2

ben1222 — Thu, 14 Oct 2010 14:16:28 +0800

Electric box 2
Electric box的第二作…
目标是通过能量的各种形式的转换将能量传输到目标点…
各种神一般的布局…

能量形式…电,光,风,流水,蒸汽,激光,无线,电池,仓鼠(这个怎么也算…)
电…有电线的地方就能传输…
光…灯泡和太阳能电池板…
风…风扇和发电机…
流水…水壶中流不完的水和发电机…
蒸汽…水壶中沸腾不完的水和…蒸汽机?
激光…激光发射器和镜面反射和激光接收器…
无线形式的电流…一次连上了无论走多远都能继续传输…
电池…充电一次就可以一直放出电…
仓鼠…看到食物或听到音乐就会跑步发电…

还有各种稀奇古怪的东西…像~~天元突破~~钻头啦…~~白井黑子~~teleporter啦…~~哈利波特与~~魔法石啦…

最后一关花了半小时才解出…
Electric box 2 Level 40

pe vol29 141~145

ben1222 — Tue, 05 Oct 2010 22:56:59 +0800

Project Euler放完假后…发现做了许多不得了的更新…
支持LaTeX这个不错…而且神奇的是…居然是以HTML形式而非图片形式显示的…真的能做到啊…这样也不会有缩放失真…真好…
然后增加了Veterans 这个等级…条件是300+…不过不会显示题数了…也就是说…没法知道是不是解决了100%的题目…
而且这个等级中的排名是按注册时间来的…
这个标准…不会改动吗?否则…比如总题数有600题的时候…貌似扯远了…
不过对防作弊而言是不错的…
最后…那个public profile页面居然不再是public的了…虽然这是引入Veterans等级后的必然结果…因为那里有写用户做了哪些题…

帖子里讨论了很多关于这次改版的事…
为什么题目提供了forum却还是很多人发到公开的博客上呢…
语言障碍和与forum中雷同的思路…很有力的观点…
反驳是…问了很多博客的主人…得到的回答大多数却是:因为别人都这么做…

自问一下…我是什么原因呢…
最初可能是因为…像在oj上做题似的…然后做题会写解题报告…然后…
不过当时貌似是因为WolframAlpha才决定做PE的…?

扯远了…
Project Euler

Problem 141 : C++,找使存在除数,商,余数,形成等比数列的完全平方数,枚举等比数列,判完全平方
Problem 142 : C++,找3个数两两求和求差都是完全平方数,枚举勾股数,注意考虑非本源勾股数
Problem 143 : C++,费马点到三边距离以及三边都是正整数的三角形,枚举内角120度的整数边三角形后配对
Problem 144 : C++,计算几何,椭圆形与直线交
Problem 145 : C++,数字反向相加后全由奇数位组成,暴力枚举

下面是详细内容:

Problem 141
Investigating progressive numbers, n, which are also square.
正整数n除以d的商是q,余数是r,如果存在d,q,r能组成等比数列那么称n是递进的
求n小于10^12中所有递进的并且是完全平方数的总和

n=d*q+r
由于r是余数,所以r肯定小于d
那么只有三种排列
r,d,q
r,q,d
q,r,d

与其枚举n,还是试试枚举等比数列吧…
公比一定是个有理数,设为R=u/v
a,b=a*R,c=a*R^2
a=v*v*k
b=v*u*k
c=u*u*k
枚举u,v,k来获得等比数列…其中u和v必须是互质的

另一个条件是n为完全平方数,看看三种排列带入后是什么样的式子
r,d,q: n=b*c+a=v*u*k*u*u*k+v*v*k=v*k*(u*u*u*k+v)
r,q,d: n=b*c+a=v*k*(u*u*u*k+v)
q,r,d: n=a*c+b=v*v*k*u*u*k+v*u*k=v*u*k*(v*u*k+1)
前两种排列得到的结果是一样的,而第三种排列的情况不可能是完全平方数(有形如x*(x+1)的完全平方数吗?)

因此只要枚举u,v,k,然后判断 v*k*(u*u*u*k+v) 是否为完全平方数就行了
限制条件:u,v,k是正整数; u>v; gcd(u,v)=1; n<10^12
判平方的时候小心溢出…
5秒搞定…答案正确

Problem 142
Perfect Square Collection
题目超短…
找最小的x+y+z使得存在整数x>y>z>0满足x+y,x-y,x+z,x-z,y+z,y-z都是完全平方数

设x-y=a^2,y-z=b^2,x-z=c^2,那么
a^2+b^2=c^2
这不就是勾股数吗…
x和y可以表示成…
y=z+b^2
x=z+c^2
其他3个完全平方数表示成…
x+y=z*2+2*b^2+a^2
x+z=z*2+a^2+b^2
y+z=z*2+b^2
那么枚举勾股数吧,对每组(a,b),枚举z,并检查x+y+z=z*3+2*b^2+a^2以得到最小值
这里写代码的时候犯了个错误…忘记考虑非本源勾股数了…
结果导致得不到结果…于是又在继续想…

设x+y=d^2,x+z=e^2,y+z=f^2
d^2=c^2+f^2
d^2=b^2+e^2
也就是说
c^2=a^2+b^2
e^2=a^2+f^2
d^2=a^2+b^2+f^2
找两组勾股数(a,b,c),(a,f,e)使得 d^2=a^2+b^2+f^2
限制条件:f>b,e>c
x=(a^2+b^2+e^2)/2
y=(e^2+b^2-a^2)/2
z=(e^2-b^2-a^2)/2

然后才发现自己一开始没考虑非本源勾股数…
改正后发现其实答案挺小的…

Problem 143
Investigating the Torricelli point of a triangle
三角形费马点:到三顶点距离和最小的点
若三角形三边和费马点到三顶点距离分别都是整数,那么称这个三角形是Torricelli三角形
求费马点到三顶点距离和小于等于120000的数的和

一开始题目看错了…以为费马点到三顶点的距离和需要是整数,但是三个距离不一定要整数…结果下面就开始另一条路上的探索了…
能用三角形三边长来计算费马点到三顶点的距离和吗?
余弦公式…
b*b=p*p+q*q-2*p*q*cos120
三边的平方用p,q,r表示…
a*a=q*q+r*r+q*r
b*b=p*p+q*q+p*q
c*c=p*p+r*r+p*r
用面积表示上面等式右边的几项…
q*r*sin120/2=q*r*sqrt(3)/4=S_TBC
q*r=S_TBC*4/sqrt(3)
p*r=S_TAB*4/sqrt(3)
p*q=S_TAC*4/sqrt(3)
q*r+p*r+p*q=S_ABC*4/sqrt(3)
等式加起来…
a*a+b*b+c*c=2*(q^2+p^2+r^2)+(qr+pq+pr)
(p+q+r)^2=p^2+q^2+r^2+2*(pq+pr+qr)
=(a*a+b*b+c*c+sqrt(3)*4*S_ABC)/2
面积用海伦公式就行了
=(a*a+b*b+c*c+sqrt(3*(a+b+c)(b+c-a)(a+c-b)(a+b-c)))/2
于是要求的”费马点到三顶点的距离和”就是…
p+q+r=sqrt((a*a+b*b+c*c+sqrt(3*(a+b+c)(b+c-a)(a+c-b)(a+b-c)))/2)
根号套根号啊…真恶心…
暴力了一些小数据出来…
a,b,c,p+q+r
43,147,152,185
49,285,296,331
57,65,73,112
73,88,95,147
95,312,343,403
97,185,208,273
111,221,280,331
127,168,205,283
147,377,437,520
152,343,387,485
193,2667,2728,2855
195,1544,1591,1729
211,1541,1560,1729

找到了这么个数列…对应上面每行第四项…

说起来三角形放大整数倍也自然是满足的
所以加个gcd(a,b,c)=1的条件
sqrt(3*(a+b+c)(b+c-a)(a+c-b)(a+b-c))
从这里找找突破口…因为里面那几项必须有3的因子才可能开根号出来
首先,如果abc都是3的倍数,那么会因为gcd(a,b,c)!=1而不会被算到
枚举了对3取余后的余数情况…得到的结论是…abc中不能有2个数都是3的倍数…
很…微妙…剪不掉多少多余的数据…

n天后…重读题目…
才发现原来p,q,r也要是正整数…
那么好办多了…
三角形TAC,TAB,TBC都是整数边三角形,还是一个角120度的…
那么跟Problem 279类似的方法生成120度角整数边三角形
然后按最短边排序
枚举看能不能把它们拼成一个三角形

不到1秒就出结果了,答案正确

Problem 144
Investigating multiple reflections of a laser beam.
椭圆形白体(white cell?我把它当作黑体的反义词可以吗?)
椭圆形的方程是4*x^2+y^2=100
开口在顶端-0.01<x <0.01的位置
一束光沿某条直线射进去,问反射多少次可以从那个洞口出来

计算几何吗…
先看看切线的方程吧…
4*x^2+y^2=100
对x求导,找斜率
8x+2*y*dy/dx=0
dy/dx=-4x/y
在(x0,y0)处的切线方程是…
4*x0*x+y0*y=100
那么椭圆与直线的交点怎么算呢…
ax+by+c=0
(x/d)^2+(y/e)^2=1
椭圆的参数方程是…
x=d*cos(t)
y=e*sin(t)
代回去就是…
a*d*cos(t)+b*e*sin(t)+c=0
变形…
b*e*sin(t)=-c-a*d*cos(t)
两边平方…
(b*e)^2*(1-cos^2(t))=c^2+2*c*a*d*cos(t)+(a*d)^2*cos^2(t)
设u=cos(t),代入上式…解出u
(b*e)^2-(b*e*u)^2=c^2+2*c*a*d*u+(a*d*u)^2
u^2*((a*d)^2+(b*e)^2)+2*c*a*d*u+c^2-(b*e)^2=0
u=(-2*c*a*d±sqrt(4*(c*a*d)^2-4*((a*d)^2+(b*e)^2)*(c^2-(b*e)^2)))/(2*((a*d)^2+(b*e)^2))
=(-c*a*d±sqrt((c*a*d)^2-((a*d)^2+(b*e)^2)*(c^2-(b*e)^2)))/((a*d)^2+(b*e)^2)
设ad=a*d,be=b*e
=(-c*ad±sqrt((c*ad)^2-(ad^2+be^2)*(c^2-be^2)))/(ad^2+be^2)
=(-c*ad±sqrt((ad*be)^2+be^4-(c*be)^2))/(ad^2+be^2)
=(-c*ad±be*sqrt(ad^2+be^2-c^2))/(ad^2+be^2)
吐…在电脑上用文本来推导真是让人想吐…
估计也不会有人看吧…也懒得弄成图片了…
得到交点后就好办了,按切线方向反射,算下一个交点

#define EPS 1e-8
#define IM(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1])
/*两点转换为解析式,ret[0]=a,ret[1]=b,ret[2]=c,ax+by+c=0*/
#define PT_LINE(p1,p2,ret) {ret[1]=p2[0]-p1[0];ret[0]=p1[1]-p2[1];ret[2]=-IM(p1,ret);}
/*
直线与椭圆的交点,返回值为交点个数
ax+by+c=0
(x/d)^2+(y/e)^2=1
x=( -c*ad*d土be*d*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
y=( -c*be*e干ad*e*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
*/
long leins(double d,double e,double*pl1,double*pl2,double*ret1,double*ret2){
 double vl[3];
 PT_LINE(pl1,pl2,vl)
 double ab=vl[0]*vl[0]*d*d+vl[1]*vl[1]*e*e,delta=ab-vl[2]*vl[2];
 if(delta<=-EPS)return 0;
 else if(delta<EPS){
 ret1[0]=-vl[2]*vl[0]*d*d/ab;
 ret1[1]=-vl[2]*vl[1]*e*e/ab;
 return 1;
 }else{
 double sq=sqrt(delta);
 ret1[0]=(-vl[2]*vl[0]*d*d+vl[1]*e*d*sq)/ab;
 ret1[1]=(-vl[2]*vl[1]*e*e-vl[0]*e*d*sq)/ab;
 ret2[0]=(-vl[2]*vl[0]*d*d-vl[1]*e*d*sq)/ab;
 ret2[1]=(-vl[2]*vl[1]*e*e+vl[0]*e*d*sq)/ab;
 return 2;
 }
}

Problem 145
How many reversible numbers are there below one-billion?
对正整数n,如果[n+reverse(n)]的每一位都是由奇数组成的,那么称n是reversible
问10^9以内有多少个这样的n
其中n和reverse(n)都不能有前置0

枚举而已吗…那么多人过…
暴力只要90秒…
更好点的方法是按位数来算

FFT 快速傅立叶变换(任意点)

ben1222 — Wed, 29 Sep 2010 17:06:14 +0800

这个月几乎都在搞保研的事情…结果没拿到资格…算了…也算经历过了…

心血来潮想搞搞FFT玩玩…
虽然在数字信号处理课上有比较详细的讲过…不过是基2的,对于不是2^N点的FFT只能补零吗…
本来是想做做图像的FFT看看…但是图像的边不是2^N点的…而且补零可能也不方便…做滤波估计还会被补的0影响…

不管怎么样…还是先看看基2的FFT怎么做吧…
稍微看了下Wiki上的…貌似还有不少别的方法的FFT…
还是先看看经典的Cooley Tukey算法吧…
离散傅立叶变换的定义式大约是这个样子的

这里用W方便表示,W(a,n)表示极角是a/n个360度的角度,这里关于角度的性质W也有

作为例子,展开个n=8的情况

为了计算右边那一列n个数,直接算要n^2级别的计算量
现在将黄色和绿色的部分分别取出来

X_e和X_o分别是原数列的偶数项和奇数项的离散傅立叶变换,递归的子结构出现了
放到原式中就变成了

这样,只要分别计算出X_e和X_o并附加2*n级别的计算量,就能完成原本要n^2级别的计算量才能完成的运算了
不过毕竟是最基础的方法…代码的运行效率还有很多优化的余地

虽然说是基2的,但是实际上只是个特例
实质上这个方法并不是对分成2份才有效
比如以n=6作为例子,既可以分成2份:

也可以分成3份:

而且依然满足递归的子结构的要求
也就是,n=a*b的情况下,可以分成a份,计算a个长度为b的子序列的离散傅立叶变换值,附加a*n的级别的运算量,就可以完成n点的离散傅立叶变换了
这样,对任意的n都可以做快速傅立叶变换了:分解质因数
可惜n是质数的时候,貌似只能用定义式了…

基2,4,8的情况下,在代码级别上还可以有很好的优化
因为这时会让中间很多步骤中有关W的复数运算简化很多,不用复数乘法,而用实部虚部做较少的加减法来完成

浮点数真是功不可没啊…用很小的精度代价让一个浮点数包含了原序列每个数的一部分信息…就像全息照相一样呢…

contest zju Monthly August 2010

ben1222 — Mon, 23 Aug 2010 15:41:06 +0800

昨天是ZOJ Monthly, August 2010
好久没做题了…N天前就看到这次月赛是东方专场…于是想着一定要去凑凑热闹…
不过最后还是只做了4个小时…最后一小时忙别的事情去了…

一看题数…居然有13题…
题号也变成了各种人物的名字了…

一上来是第阿求题…禁语吗……glob和regex的转换…看了下样例感觉不难…但是害怕有一堆要考虑的东西…

然后是第⑨题…算术教室…n个人围成一圈随机选取m个人,其中有9人是连着坐的概率是多少呢?…第一反应就是感觉是组合数方面的于是无视中…

然后看到已经有人过题了…第妖梦题…为大胃王幽幽子准备食物…第i天能准备bi份食物,但是要吃掉ai份食物…希望食物越新鲜越好…求各天准备多少食物…果然很水呢…从后往前扫一遍就可以了…
本想就看看题目娱乐下…结果一不小心手痒了…
于是敲了下…过了…

接着发现又一题有人过了…第四季题…黑白分明呢…用了Bad Apple影绘的图…
把24位点阵图转换成灰度然后再根据公式二值化…也很水…看到那句”All divisions are truncated divisions.”后就想…会有人栽在这里吧…后来真发现有个id是”一句All divisions are truncated divisions. 让我wa了10次，杯具了整个下午！！！”
但是敲了一下…在本地测试的时候用scanf(“#%2lx%2lx%2lx”,&r,&g,&b);来读入,结果本地测试死循环了…难道#也是特殊字符?于是改成了scanf(“%#%2lx%2lx%2lx”,&r,&g,&b);在本地测试通过了…可是交上去OLE了…难道编译器不同吗…
于是很郁闷的改成了scanf(“%s”,str);sscanf(str+1,”%2lX%2lX%2lX”,&r,&g,&b);
题目描述里面没说输出的9或者空格之间需要空格隔开…但是sample output里面却隔开了…于是改成没有空格隔开结果wa了…
最后再改回来有空格隔开的才AC…

接着看到第灵梦题有人过了…塞钱箱…说是第i天请人赛钱可以得到si元,但是如果这么做了,下一次请人赛钱的时间必须在xi天后,yi天前…求如果第一天请人赛钱的话,最多能得到多少钱…
稍微想了下,一个简单的dp而已…dp[i]表示第i天请人赛钱的话,从第i天开始最多能得到多少钱…转移时只要找xi天后,yi天前这段区间的dp值最大的即可…
然后一看规模…5万天…囧了…不过求区间最大值而已…解题报告虽然说用线段树…可当时我想到的却是rmq…
但是rmq是一次性的预处理的…于是就YY着改造了一下,依然沿袭rmq那种思路:分成logN层,每层有N个数,第i层的第j个数表示在[j,j+2^i-1]的范围内最大值是多少
然后从后往前边dp边insert,对于第i天来说,这天之后的每一天的dp值和rmq数组值都是已经得到了的,于是就能得到第i天的dp值,然后将这个值插入进rmq数组,更新第i个数的那logN层数字,总复杂度O(nlogn)
很快写好了…提交却下标越界…然后发现数组开小了- -|||…改后AC
在看解题报告的时候发现评论里有说Sparse Table什么的…原来我那方法叫这个名字啊…

接下来看到第辉夜题有人过…妹红要埋伏辉夜吗…题目是求一张图上从起点到终点之间的必经之路…也就是桥?…模板题吗…可惜我只有割点的模板…于是忽略了这题…

接下来看到第芙兰题有人过了…永夜抄吗……看到迷途竹林以为又是图论有关的…结果却不是这样…
一共有n个房间,第i个房间有xi个ITEM,每个是ai点,有yi个TIME,每个是bi点,每次去一个房间把所有这个房间的ITEM和TIME收完回到入口,直到所有房间的ITEM和TIME都收完…
拿一个点数是a的ITEM的时候可以得到IV分,并且TV增加a
拿一个点数是b的TIME的时候可以得到TV分,并且IV增加b
一开始IV和TV值为0,求最多能得到多少分
因为进一个房间后,该房间所有的ITEM和TIME都要收完,所以先考虑在一个房间内的情况,该按什么顺序收最好
小数据手动演算了一下,发现分数会增加IV*xi+TV*yi+u*ai+v*bi,并且u+v=xi*yi,并且对任何正整数的u,v搭配都总有一种顺序可以满足
于是得到了在第i个房间内最多分数会增加IV*xi+TV*yi+max{ai,bi}*xi*yi的结论,可见ai,bi和IV,TV是分离的,因此外面用一层状态压缩dp计算2^15个状态就可以了
dp[state]表示state的位为1的对应的房间已经去过的情况下,所得到的最大分数
并且这个情况下的IV和TV值是固定的,轻轻松松的上吧…

又看到第因幡题有人过…腹黑兔啊…看似公平的游戏却把规则定的自己总有必胜策略…想到东方夜神雪里面那个取硬币小游戏…轮流取1~3枚,取走最后一枚的人输…可是起始的硬币数量并不是随机的…而是4k+1…太囧了…玩了好几盘才发现原来这游戏是必输的…
回到这里…这次是类似孔明棋的小游戏…每次可以选一个棋子往4个方向之一跳,但是必须跳过别的棋子并且最终落在某个空格上,中间可以跳过连续多个棋子,被跳过的棋子都会被拿走…
给定初始棋局,问让铃仙和因幡谁先走可以使得因幡有必胜策略…
一看规模才5*5…想想2^25也不是什么大数字…加上有一句”It’s guaranteed that no more than 61.8% of them are ‘e’s.”…感觉状态会很少…于是直接暴力记录博弈状态…
第一次交忘记处理了选取的棋子是一排棋子中间的某个棋子的情况…改后AC

然后…第魔里沙题有人过…魔炮啊…是个计算几何…求往什么方向射可以打中最多点,魔炮的范围是个顶点为原点的抛物线
抛物线啊…想到每个点有个可被命中的角度范围,也就是魔炮的方向在这个范围内,这个点可以被打中,但是对不同距离的点而言,这个范围的大小是不一样的…
计算了下,二次系数是a的抛物线和半径为r的圆的交点

y坐标是: y=(sqrt(1+4*a*a*r*r)-1)/(2*a)

设上图中P点对应的极角为alpha,那么当魔炮中轴对着[alpha-theta, alpha+theta]区间中任意的极角时,P点就能被打中

这个theta值为: theta=atan(x/y) =atan(sqrt(y/a)/y)=atan(sqrt(1/y*a))

于是问题转换成了,找一个极角,使得它被最多的区间覆盖…而数据规模有3万个点…
脑中第一个闪现的是离散化+线段树…显然太麻烦了…
然后想起以前某些题…比如矩形求并之类的…貌似只要看左右边界判断出入就行…
于是想到了用这个方法就好,将区间变成端点,左端点标记为入,右端点标记为出,然后将所有端点按大小(这里是极角)排序
然后从前往后扫描,遇到标记为入的端点就当前覆盖区间数加一,遇到标记为出的端点就当前覆盖区间数减一,扫描一遍下来,取当前覆盖区间数的最大值…就是所需的答案daze
另外处理下跨边界相邻的情况,就是当区间在[0,2pi]的两边时,就多插入个小于0的那个入端点和大于2pi的那个出端点

然后…有别的事情…于是就看看题目然后不继续做了…
第极题…姬海棠也有啊…貌似是找最安全点…话说游戏里寅丸星的香蕉激光实在是太恶心了…但是这题目里面是射线…而且是曼哈顿距离?…直接傻眼了…点和线的曼哈顿距离?那是什么?怎么求?

第图书题…不动的大图书馆…题目是一共m种元素卡,一共有n种符文,每种元素卡上等概率出现这n种符文中的一种,求m种元素卡各取一张的话,有多少概率可以使得这些卡上的符文至少有l个是相同的…输出成既约分数…
又是组合数类的题目吗…

第咲夜题…时间的操作…哦不…这题是求表达式结果是某个类型的可能性个数?
题目很长呢…

第幽幽子题…赏花会啊…真好呢…
又是做食物…这次是第i天幽幽子会吃掉ai份食物…妖梦每天可以做L份食物,或者花一天时间锻炼自己使得L值增加1…求N天过后剩下的食物最多还有多少份…

这里或者这里是官方解题报告吧…

下面是我灵梦,芙兰,魔里沙题的代码…

第灵梦题的代码:

#include<stdio.h>
/*
ZOJ monthly 201008
红白
*/
#define N 50010
long dp[N];
//dp[i]=max_money get if ask for money at day i
long rmq[18][N],bi[19];
long n;
long money[N],begin[N],end[N];
long query(long pa,long pb){
 long j,d;
 if(pa>pb){j=pa;pa=pb;pb=j;}
 if(pa>=n||pa<0)return 0;
 if(pb>=n)pb=n-1;
 d=pb-pa+1;
 for(j=0;bi[j]<=d;j++);
 j--;
 pa=rmq[j][pa];
 pb=rmq[j][pb-bi[j]+1];
 if(pa<pb)pa=pb;
 return pa;
}
void insert(long p,long val){
 long j;
 for(j=0;j<18;j++){
  rmq[j][p]=val;
  long next=bi[j]+p;
  if(next<n){
   if(rmq[j][next]>val)val=rmq[j][next];
  }
 }
}
int main(){
 long i;
 for(i=0;i<19;i++)bi[i]=((long)1)<<i;
 while(scanf("%ld",&n)!=EOF){
  for(i=0;i<n;i++)scanf("%ld%ld%ld",money+i,begin+i,end+i);
  long j;
  for(j=0;j<18;j++)for(i=0;i<n;i++)rmq[j][i]=-1;
  for(i=n-1;i>=0;i--){
   dp[i]=money[i]+query(i+begin[i],i+end[i]-1);
   insert(i,dp[i]);
  }
  printf("%ld\n",dp[0]);
 }
 return 0;
}

第芙兰题的代码:

#include<stdio.h>
#include<memory.h>
/*
ZOJ monthly 201008
二小姐

in part_i:
 point+=IV*xi+TV*yi+max{ai,bi}*xi*yi
 IV+=yi*bi
 TV+=xi*ai
dp[state(1<<n)]=max_point when collected in parts[state.bits]
*/
#define N 32768
long bi[17];
long dp[N],s_iv[N],s_tv[N];
long n;
long xi[15],yi[15],a[15],b[15];
long add_point[15];
long fun(long state){
 if(dp[state]!=-1)return dp[state];
 long i;
 long best=0;
 for(i=0;i<n;i++)if(state&bi[i]){
  long prev=state^bi[i];
  long iv=s_iv[prev],tv=s_tv[prev];
  long val=fun(prev)+iv*xi[i]+tv*yi[i]+add_point[i];
  if(val>best)best=val;
 }
 return dp[state]=best;
}
int main(){
 long i;
 for(i=0;i<17;i++)bi[i]=((long)1)<<i;
 while(scanf("%ld",&n)!=EOF){
  for(i=0;i<n;i++){
   scanf("%ld%ld%ld%ld",a+i,b+i,xi+i,yi+i);
   if(a[i]>b[i])add_point[i]=a[i]*xi[i]*yi[i];
   else add_point[i]=b[i]*xi[i]*yi[i];
  }
  memset(dp,-1,sizeof(dp));
  long j;
  for(i=0;i<bi[n];i++){
   long iv=0,tv=0;
   for(j=0;j<n;j++)if(i&bi[j]){
    iv+=yi[j]*b[j];
    tv+=xi[j]*a[j];
   }
   s_iv[i]=iv;
   s_tv[i]=tv;
  }
  printf("%ld\n",fun(bi[n]-1));
 }
 return 0;
}

第魔里沙题的代码:

#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
/*
ZOJ monthly 201008
黑白

y=a*x^2
x*x+y*y=r*r
y=a*(r^2-y^2)
y=a*r*r-a*y*y
a*y^2+y-a*r*r=0
y=(sqrt(1+4*a*a*r*r)-1)/(2*a)
theta=atan(x/y) =atan(sqrt(y/a)/y)=atan(sqrt(1/y*a))
*/
typedef struct{
 double rad;
 char out;
}point;
bool operator<(const point&a,const point&b){
 return a.rad<b.rad||(a.rad==b.rad&&a.out<b.out);
}
#define N 30010
#define PI 3.14159265358979
point seg[N*3];
double xi[N],yi[N];
long nseg;
#define ENSEG(r,o) {seg[nseg].rad=r;seg[nseg].out=o;nseg++;}
int main(){
 long n;
 double a;
 while(scanf("%ld%lf",&n,&a)!=EOF){
  long i;
  for(i=0;i<n;i++)scanf("%lf",xi+i);
  for(i=0;i<n;i++)scanf("%lf",yi+i);
  long addition=0;
  nseg=0;
  for(i=0;i<n;i++){
   double r2=xi[i]*xi[i]+yi[i]*yi[i];
   double dy=(sqrt(1.0+4*a*a*r2)-1)/(2*a);
   if(dy<=0){addition++;continue;}
   double theta=atan(sqrt(1.0/(dy*a)));
   double t0=atan2(yi[i],xi[i]);
   if(t0<0)t0+=PI*2;

   double in=t0-theta,out=t0+theta;
   ENSEG(in,0)
   ENSEG(out,1)
   if(in<0){
    ENSEG(in+PI*2,0)
    ENSEG(out+PI*2,1)
   }
   if(out>=PI*2){
    ENSEG(in-PI*2,0)
    ENSEG(out-PI*2,1)
   }
  }
  sort(seg,seg+nseg);
  long ans=0,now=0;
  for(i=0;i<nseg;i++){
   if(seg[i].out)now--;
   else now++;
   if(now>ans)ans=now;
  }
  printf("%ld daze\n",ans+addition);
 }
 return 0;
}

objective-c 的一些东西

ben1222 — Fri, 13 Aug 2010 23:03:48 +0800

最近在看Objective-c
虽说是c语言扩展出来的…可是这里面的对象甚至类本身却是动态的…
查了下貌似这叫反射?
倒是和不少脚本型的语言很像…

了解下语法之类的可以看看这里(这篇东西同时有中文版),Wiki上还有些和C++的对比

一点点看下来,一开始还没什么,感觉只是语法变了下
嗯…#include变成#import

创建类用@interface,还可以继承,不过只能继承一个父类
objective-c的方法名很有特点…并不是完整名称后面跟参数列表…而是两者交织起来…
虽然这样打字可能要多些(没有带提示功能的编辑器的话,很吃力…),但是阅读起来更加明晰,不会出现忘记第几个参数是什么作用的情况了

实现类在@implementation块内写
这里就有些舒服了,方法名直接复制到对应的implementation块中就可以开始写方法的内容了
用C++的时候,要把成员函数的声明复制过来,然后在返回类型和函数名之间加上”类名::”的东西,感觉很麻烦呢…

this变成了self,要调用父类的方法而自身却重载过了这个方法时,可以对super调用方法
用-的是对象的method,用+的是类的method,后者相当于C++的静态成员函数

id类型是万能的对象,就像void*,对id对象可以调用任何方法,只要它的本体有这个方法就可以,没有的话…就会产生异常
看到后面感觉越来越玄乎…难道一个对象里存了一个方法列表,并且方法名和参数类型都有吗?…
事实上就是这样没错…一个对象指针所指的是个这样的结构体(摘自objc.h)

typedef struct objc_object {
  struct objc_class*  class_pointer;
} *id;

实际对象在这个基础上后面还会跟上它的成员变量…
而这个objc_class呢,则是这样的结构(摘自objc.h)

struct objc_class {
  MetaClass           class_pointer;          /* Pointer to the class's
                                                meta class. */
  struct objc_class*  super_class;            /* Pointer to the super
                                                class. NULL for class
                                                Object. */
  const char*         name;                   /* Name of the class. */
  long                version;                /* Unknown. */
  unsigned long       info;                   /* Bit mask.  See class masks
                                                defined above. */
  long                instance_size;          /* Size in bytes of the class.
                                                The sum of the class
						definition and all super
						class definitions. */
  struct objc_ivar_list* ivars;               /* Pointer to a structure that
                                                describes the instance
                                                variables in the class
                                                definition.  NULL indicates
                                                no instance variables.  Does
                                                not include super class
                                                variables. */
  struct objc_method_list*  methods;          /* Linked list of instance
                                                methods defined for the
                                                class. */
  struct sarray *    dtable;                  /* Pointer to instance
					         method dispatch table. */
  struct objc_class* subclass_list;           /* Subclasses */
  struct objc_class* sibling_class;

  struct objc_protocol_list *protocols;	      /* Protocols conformed to */
  void* gc_object_type;
};

其中值得注意的是objc_method_list,里面果然存在一个方法列表,只要在列表里查询匹配的方法就可以了,找不到还能在父类中找
再稍微深入一下,看看这个objc_method_list是什么东西(摘自objc-api.h)

typedef struct objc_method_list {
  struct objc_method_list*  method_next;    /* This variable is used to link
                                               a method list to another.  It
                                               is a singly linked list. */
  int            method_count;              /* Number of methods defined in
                                               this structure. */
  Method method_list[1];                    /* Variable length
                                               structure. */
} MethodList, *MethodList_t;

这是一个…链表…而且似乎是块状链表?
继续,找到Method的定义是(摘自objc-api.h)

typedef struct objc_method {
  SEL         method_name;                  /* This variable is the method's
                                               name.  It is a char*.
                                               The unique integer passed to
                                               objc_msg_send is a char* too.
                                               It is compared against
                                               method_name using strcmp. */
  const char* method_types;                 /* Description of the method's
                                               parameter list.  Useful for
                                               debuggers. */
  IMP         method_imp;                   /* Address of the method in the
                                               executable. */
} Method, *Method_t;

可以看到,一个方法包含一个SEL,一个参数表,一个IMP
SEL就是用@selector()可以得到的东西,上面虽然说包含方法名,但是我实际测试下来似乎不是这样…?版本不同吗?
IMP是真正可执行的函数的地址

于是像什么respondsToSelector:(判断一个对象有没有某个方法),performSelector:(执行某个方法),poseAsClass:(一个类完全由其某个子类扮演),这种神奇的函数也就不那么神秘了…

在objective-c中,类自身也是运行时存在的一个对象,想获取这个类的信息(比如类名)或修改这个类,应该是可以做到的
要新建一个类的实例,需要用类的alloc方法
也就是说,与其说是类,不如说是一个工厂对象…生产实例,以及记录这个实例所支持的方法列表之类的信息

做了个实验,建了个不继承其他类的类BObject,继承自BObject的类Rectangle
建了BObject,Rectangle,NSObject的实例,然后查看它们的指针以及meta类和父类的指针,画了下面的图

从方法列表的方法数来看,似乎左侧的只包含了”-”方法,其meta类包含了”+”方法
上面实验是在ubuntu下的GNUStep环境中做的…

稍微研究下机理,可以帮助了解很多本来不怎么好理解的特性,函数,用法,等

这里还有几篇较深入的文章

pe vol28 136~140

ben1222 — Tue, 03 Aug 2010 22:14:51 +0800

原来project euler也是会放假的…出到299题后停2个月…说会有新的等级出现…

还是慢慢啃老题吧…

Project Euler

Problem 136 : C++,类似Problem 135
Problem 137 : C++,att,考察菲波那契数列为系数的无穷级数收敛到正整数的情况,小数字搜到大数字…
Problem 138 : C++,att,考察边是整数的等腰三角形高与底边大小相差1的情况,还是小数字搜到大数字…
Problem 139 : C++,直角三角形拼正方形,观察一下实际上所需满足的条件,然后枚举本源勾股数处理一下就行了…
Problem 140 : C++,类似Problem 137,使用广义菲波那契数列,小数字找规律得到递推式…

下面是详细内容:

Problem 136
Discover when the equation x² − y² − z² = n has a unique solution.
正整数x,y,z是等差数列,对给定的正整数n,判断x^2-y^2-z^2=n是否只有唯一解
求n<50,000,000中上式只有唯一解的n的个数

和Problem 135很像啊…
用相同的方法试试…
只不过从求f(n)=10的个数变成了f(n)=1的个数
6秒出结果…

Problem 137
Determining the value of infinite polynomial series for which the coefficients are Fibonacci numbers.
无穷级数AF(x)=x*F1+x²*F2+x³*F3+…
其中Fk是第k个菲波那契数,F1=1,F2=1
求使AF(x)为正整数的第15个有理数x所对应的AF(x)的值

先看看级数收敛成什么样吧…

AF(x)/x - AF(x) = F1+x*(F2-F1)+x²*(F3-F2)+...
                = F1+x*F0+x²*F1+x³*F2+...
因为F1=1,F0=0
AF(x)/x - AF(x) - 1 = x²*F1+x³*F2+...
                    = x*AF(x)
AF(x)*(1/x-1-x)=1
AF(x)=1/(1/x-1-x)=x/(1-x-x²)

原来是这么一个式子…
那么现在需要x是有理数,AF(x)是正整数…

设AF(x)=y
y*(1-x-x²)=x
x²*y+x(y+1)-y=0
x=(-(y+1)+sqrt((y+1)²+4*y*y))/(2*y)

使x为有理数就必须让根号内的数是完全平方数
(y+1)²+4*y*y=(y+1)²+(2*y)²
啊…好像勾股数耶…
不过按本源勾股数来的话结果漏掉一些有效解…而且数字太大…吃不消…

于是换个方法,先枚举测试较小值看看有什么方法…

1:2
2:15
3:104
4:714
5:4895
6:33552
7:229970
8:1576239
9:10803704
10:74049690

拿去att搜索一下发现有个叫做golden rectangle的东西…
啊…难道就是这个数列的偶数项吗…?
这种数的公式是:F(n)*F(n+1)

然后又发现了这个数列
答案居然直接在上面…

Problem 138
Investigating isosceles triangle for which the height and base length differ by one.
等腰三角形,底边为b,腰为L,高为h
求最小的12个等腰三角形的L的和,满足h=b土1
其中b,L,h为正整数

总之先搞来勾股数然后一一判断看看…
两直角边如果满足一边与另一边的两倍的差是1的话就行了
似乎底边,也就是b,必须是偶数?

发现这数的增长速度真快…
拿L的前几项去att搜了一下…发现是这个数列

理由呢?
好像这种数有个特性…5*L²-1是完全平方数…
验证下是不是满足题目中的三角形…

设a=b/2
那么由等腰三角形的性质可得:
L²=a²+(2*a土1)²
=5*a²土4*a+1
5*L²-1=25*a²土20*a+4=(5*a土2)²

啊…果然如此…

Problem 139
Finding Pythagorean triangles which allow the square on the hypotenuse to be tiled.
用四个直角三角形拼成一个正方形,外加中间一个正方形的洞
用中间的洞的大小的小正方形能填充大正方形的区域的话,就算成功
问周长在100,000,000以内的直角三角形有几个能成功

(a,b,c)用c边拼大正方形,中间的洞的正方形边长是(b-a)
只要c能被(b-a)整除就能成功
生成勾股数然后判断就行了
只要注意到本源勾股数能成功的话,它的任何倍数都会成功
而本源勾股数不能成功的话,它的任何倍数都不会成功

Problem 140
Investigating the value of infinite polynomial series for which the coefficients are a linear second order recurrence relation.
类似Problem 137
无穷级数AG(x)=x*G1+x²*G2+x³*G3+…
其中Gk是一个广义菲波那契数列的第k个数,G1=1,G2=4
求使AG(x)为正整数的前30个有理数x所对应的AG(x)的值的和

老规矩…看看级数如何收敛

AG(x)/x - AG(x) = G1+x*(G2-G1)+x²*(G3-G2)+...
                = G1+x*G0+x²*G1+x³*G2+...
因为G1=1,G0=3
AG(x)/x - AG(x) - 1 = x*3+x²*G1+x³*G2+...
                    = x*3+x*AG(x)
AG(x)*(1/x-1-x)=1+x*3
AG(x)=(1+x*3)/(1/x-1-x)=(x+3*x²)/(1-x-x²)

那么现在需要x是有理数,AG(x)是正整数…

令AG(x)=y
y*(1-x-x²)=x+3*x²
x²*(3+y)+x*(1+y)-y=0
x=(-(1+y)+sqrt((1+y)²+4*y*(3+y)))/(2*(3+y))
根号内要是完全平方数
(1+y)²+4*y*(3+y)=1+2*y+y²+12*y+4*y²
                 =5*y²+14*y+1

枚举测试较小数的话…

这回在att上搜不到了,但是因为Problem 137最后是菲波那契数列相乘的形式,所以这回试着对这些数做做差找找规律
结果很快就发现了这样的规律
y3=y2+G5+F3
y4=y3+G6-F3
y5=y4+G9+F7
y6=y5+G10-F7
y7=y6+G13+F11
y8=y7+G14-F11
..
其中Fk是第k个菲波那契数,F1=1,F2=1

剩下的就是递推了,答案正确