Project Euler 291
@see http://projecteuler.net/index.php?section=problems&id=291
/**
* Copyright uwi ( http://wonderfl.net/user/uwi )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/zqpC
*/
package {
import flash.display.Sprite;
import flash.text.TextField;
import flash.utils.getTimer;
// @see http://projecteuler.net/index.php?section=problems&id=291
public class Euler291 extends Sprite {
private var _tf : TextField;
// p=(x^4-y^4)/(x^3+y^3) (x,yは正の整数)の形で表される
// 5*10^15未満の素数pの個数はいくつか
public function Euler291() {
_tf = new TextField();
_tf.width = 465;
_tf.height = 465;
addChild(_tf);
var s : int = getTimer();
tr(solve(5E12));
// tr(solve(5E15));
var g : int = getTimer();
tr((g - s) + " ms");
}
// まず明らかに、x-y>0, xy>0
// (x^4-y^4)/(x^3+y^3)=(x-y)(x^2+y^2)/(x^2-xy+y^2)
// =AB/Cとすると、
// 0<A<=A^2<C<Bより、
// AB/Cが素数になるためには、CがAで割り切れなければならない。
// x-y|x^2-xy+y^2 -> x-y|xy -> xy=K(x-y) (Kは自然数)
// ここで、AB/C=(x-y)(x-y+2K)/(x-y+K) (=素数)
// となる。x-y < x-y+K < x-y+2Kより、
// x-y+Kがx-yで割り切れなければならない。
// x-y|x-y+K -> x-y|K -> K=L(x-y) (Lは自然数)
// とすると、AB/C=(x-y)(1+2L)/(1+L) (=素数)
// となる。1+2Lと1+Lは互いに素なので、
// AB/C=2L+1, x-y=1+Lでなければならない。
// このとき、xy=L(x-y)^2=L(1+L)^2なので、
// 2次方程式の解と係数の関係より、
// -x,y = (-(1+L)±√((1+L)^2+4L(1+L)^2))/2
// = (1+L)(-1±√(1+4L))/2
// となるので、1+4Lが平方数であればx,yが存在する。
// 4L+1は奇数なので、自然数tを用いて4L+1=(2t+1)^2とかける。
// よって、L=t(t+1)となるので、AB/C=2t(t+1)+1となる。
// あとは、自然数tに対して、p(t)=2t(t+1)+1が素数かどうか判定すればいいのだが、
// p(t)<=5*10^15のとき、t<=5*10^7なので、地道にやっていくと
// かなりの時間がかかる。
//
// 補題1 v|p(t) (vは素数)のとき、整数kに対し必ず
// v|p(t+kv), v|p(-t-1+kv)
// が成り立ち、それ以外のpはvの倍数ではない。
// 証明 p(w)-p(t)=2(w-t)(w+t+1)より、これがvの倍数になる
// ことは、w-tがvの倍数またはw+t+1がvの倍数であることと同値なのでOK.
//
// 補題2 p(u)の√p(u)以下の任意の素因数について、
// あるt<uが存在してp(t)の素因数にもなっている。
// 証明 補題2が成り立たないと仮定する。
// p(u)の√p(u)以下の素因数をdとすると、
// d < √(2u(u+1)+1) < √2*u+1
// 補題1より、d|p(u-d)とならないように u-d<=0
// また、d|p(d-u-1)とならないように、d-u-1<=0 or d-u-1>=u
// 以上の3つの条件を合わせると、d=u or u+1となるが、
// どちらの場合もp(u)を割り切らないので不適。
//
// 補題2より、p(u)を、p(u)より前のpの素因数で割れるだけ割ってしまうと、
// 1もしくは素数になるため、これを利用して、素数判定をしなくても
// 逐次求めていける。
private function solve(N : Number) : Number
{
var sup : uint = (Math.sqrt(2*N-1)-1)/2;
var map : Vector.<Number> = new Vector.<Number>(sup+1);
var t : uint;
for(t = 1;t <= sup;t++){
map[t] = 2*t*t+2*t+1;
}
var ct : uint = 0;
var i : uint;
var p : uint;
for(t = 1;t <= sup;t++){
var n : Number = map[t];
if(n == 1)continue;
if(n == 2*t*t+2*t+1){
ct++;
}
for(i = t;i <= sup;i += n){
while(map[i] % n == 0)map[i] /= n;
}
for(i = n-1-t;i <= sup;i += n){
while(map[i] % n == 0)map[i] /= n;
}
}
return ct;
}
private function sieveEratosthenes(n : int) : Array
{
var nn : uint = ((n / 2 - 1) >>> 5) + 1;
var ar : Vector.<uint> = new Vector.<uint>(nn);
var i : uint, j : uint;
for(i = 0;i < nn - 1;i++)ar[i] = 0xffffffff;
ar[nn - 1] = (1 << ((n / 2 - 1) & 31)) - 1;
var sq : uint = (Math.sqrt(n) - 3) >>> 1;
for(var p : uint = 0;p <= sq;p++){
if(ar[p >>> 5] & (1 << (p & 31))){
var m : uint = (p << 1) + 3;
var m2 : uint = m << 1;
for(var mm : uint = m * m;mm <= n;mm += m2){
var ind : uint = (mm - 3) >>> 1;
ar[ind >>> 5] &= ~(1 << (ind & 31));
}
}
}
var ret : Array = [2];
for(i = 0;i < nn;i++){
for(j = 0;j <= 31;j++){
if(ar[i] & (1 << j))ret.push((((i << 5) | j) << 1) + 3);
}
}
return ret;
}
private function factor(n : Number) : Object
{
var ret : Object = {};
var sq : int = Math.sqrt(n);
for(var i : int = 2;i <= sq;i++){
var ct : int = 0;
while(n % i == 0){
n /= i;
ct++;
}
if(ct > 0){
ret[i] = ct;
sq = Math.sqrt(n);
}
}
if(n != 1){
ret[n] = 1;
}
return ret;
}
private function isPrime(n : Number, t : uint = 3) : Boolean
{
if(n <= 1)return false;
if(n % 2 == 0)return n == 2;
// if(_cache[n])return _cache[n] == 1;
var sq : int = Math.sqrt(n);
// t += (t-1) % 2;
for(var i : int = 3;i <= sq;i+=2){
if(n % i == 0){
// _cache[n] = -1;
return false;
}
}
// _cache[n] = 1;
return true;
}
private function tr(...o : Array) : void
{
_tf.appendText(o + "\n");
_tf.scrollV = _tf.maxScrollV;
}
}
}