プロジェクトオイラー問94 - sinapusu2002 @ ウィキ

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2094
Problem 94 「擬正三角形」 †
一辺の長さが整数の正三角形は面積が整数にならないことを示すのは簡単である. しかし, 5-5-6の辺を持つ殆ど正三角形に近い擬正三角形 (almost equilateral triangle) は面積が12で整数である.

以降, 二等辺三角形で, 3つめの辺の長さが他と1つしか違わないもの (5-5-6, 5-5-4等) を, 擬正三角形と呼ぶ.

さて, 周囲の長さが1,000,000,000以下の面積が整数になる擬正三角形を考え, その周囲の長さの総和を求めよ.

解法
ヘロンの公式の式変形から。全探索で求めてみました。
４０秒もかかっています。
遅いです。
なんか正当なコードは行列演算ででるみたいです。
原理がわからないのでこういう時無学な自分を痛感します。

#include<stdio.h>
#include<iostream>
#include <time.h> 

bool ok16A(__int64 a){
	__int64 a1=3*a+1;
	int i16=16;
	while(a1%2==0){
		a1/=2;
		i16/=2;
	}
	a1=a+1;
	while(a1%2==0){
		a1/=2;
		i16/=4;
	}
	a1=a-1;
	while(a1%2==0){
		a1/=2;
		i16/=2;
	}
	return i16<2;
}
bool ok16B(__int64 a){
	__int64 a1=3*a-1;
	int i16=16;
	while(a1%2==0){
		a1/=2;
		i16/=2;
	}
	a1=a-1;
	while(a1%2==0){
		a1/=2;
		i16/=4;
	}
	a1=a+1;
	while(a1%2==0){
		a1/=2;
 		i16/=2;
	}
	return i16<2;
}
int main(){
	clock_t start = clock();    // スタート時間
	__int64 m=2,n=1,ans=0,n3;
	//a a a+1を解く
	while(1){
		__int64 m2=m*m;
		__int64 n2=4+3*n*n;
		__int64 n3=m+1;
		if(m2==n2){
			ans+=ok16A(n3)*(n3+1);
			m+=3;
		}else if(m2>n2){
			n++;
		}else{
			m+=3;
		}
		if(n3+1>=10000*100000)break;
		if(n%10000000==0)std::cout<<n<<" "<<ans<<"\n";
	}
	m=4,n=1;
	std::cout<<",,,,,,,,\n";
	while(1){
		__int64 m2=m*m;
 		__int64 n2=4+3*n*n;
		__int64 n3=m-1;
		if(m2==n2){
			ans+=ok16B(n3)*(n3+1);
			m+=3;
		}else if(m2>n2){
			n++;
		}else{
			m+=3;
		}
 		
		if(n3-1>=10000*100000)break;
		if(n%10000000==0)std::cout<<n<<" "<<ans<<"\n";
	}
	std::cout<<"ans="<<ans<<"\n";
	 clock_t end = clock();     // 終了時間
       std::cout << "duration = " << (double)(end - start) / CLOCKS_PER_SEC << "sec.\n";

}