2.	実数に関するプログラム

2.1	機械イプシロン

	double machineEP(){
		double E = 1, ED;
		do{ ED = E; E = E / 2;} while((E + 1) > 1);
		E = 1 + ED; return E - 1;
	}

2.2	自前のPower

	double myPower(double X,int N){
		double P=1.0;
		for(int abN = abs(N); abN > 0; abN >>=1, X *= X) if(abN & 1) P *= X;
		if(N < 0) return 1/P; else return P;
	}

2.3	実数を連分数に変換する
(1)
	int continuedFraction(float X,int N,int A[]){
		double B = floor(X); A[0]=(int)B;
		if(fabs(B - X) < 0.00000001) return 1;
	    	for(int i=1;i<N;i++){　X = 1 / (X - B); B = floor(X); A[i]=(int) B;	}
		return 20;
	}
(2)
	double contFracToDouble(int N, int A[]){
		double X = 0; for(int i = N - 1; i>0;i--) X = 1/(A[i]+X);
		return X + A[0];
	}

2.4	自前のSin関数
	double mySin(double X){
  		double bS, S, H, DH, Pi, HPi, XX, A; int K, i;
  		H  = 3.23846264338328 / 2; DH = 1E-15  ; 
  		Pi = 3.14159265358979    ; HPi = Pi / 2;
  		XX = fabs(X); K = (int) floor(XX / HPi); /* π/2の整数倍を求め差し引く*/
  		XX = XX - ((double) K) * HPi -((double) K) * H * DH; /* 桁落ち対策     */
  		A  = -XX * XX; i = 1;                    /* π/2の奇数倍のとき　       */
  		if ( K & 1){ i = 0; XX = 1; }            /* cosを求めるものとする      */
  		if ( K < 0) K --;                        /* K<0のときのFix値調整      */
  		if ((K / 2) & 1) XX = -XX;               /* k / 2が奇数のとき符号逆転  */
  		S = XX;                      /* マクローリン級数の初期値をXXとする     */
  		do{                          /* マクローリン展開の収束計算              */
    			bS = S; i += 2; XX *= A / ((double)((i - 1) * i)); S += XX;
  		} while (S != bS);           /* 加算結果が変わらないとき完了            */
 		if( X < 0) S= -S ;            /*  元の値が負のとき符号逆転。             */
 		return S;
	}

2.5	自前のExp関数
	double myExp(double X){
 		int K, N; double E, AX, XX, bE, Alf, P, PP, A1, A2;
 		XX = fabs(X); K=(int)XX; XX -= (double)K; /* 整数分を差し引く */
 		E = 1; N = 1; AX = 1;                     /* 初期設定 */
 		do{                                       /* Expのマクローリン展開の収束計算*/
   			bE = E; AX *= XX / (double)N; E += AX; N++;
 		} while (E != bE);                       /* 加算結果が変わらないとき完了*/
 		Alf = 0.523536028747135; P  = Alf / (1 - Alf); /* 補正用設定*/
 		for(PP = 1 + P, A1  = 2.71828182845904, A2 = 2.71828182845905;
     			K > 0; K >>= 1, A1 *= A1, A2 *= A2)  /* 回数を2で除し乗ずる値を2乗する。*/
     		if( K & 1) E *= (A1 + P * A2) / PP;  /* 2進数該当桁が1のとき乗算（補正） */
		if( X < 0) return 1 / E;                /* 元の値が負のとき逆数にする。*/
 		return E;   
	}

2.6	連分数による関数計算
(1)
	double myLog(double X){
 		int N = 40;double log2, Flog2;     /* Flog2は桁落ち対策用*/
 		int i;double XX, A, L;
 		log2 = 0.693147180559945; Flog2 = 3.09417232E-16;
 		XX = fabs(X); if(XX < 0.0000000000001) return -1E+308;
 		L = 0;
 		while(XX >= 2){                    /* log2分を加算           */
   			XX /= 2; L = L + log2 + Flog2;   /* ここで桁落ちはないが       */
 		}                                  /* 減算のときと式を揃えるため */
 		while(XX < 1){                     /* XX < 1 のとき差し引く*/
   			XX *= 2; L = L - log2 - Flog2;   /* 桁落ち対策 */
 		}  /*  上記で L -= (log2 + Flog2)では桁落ち対策にならないことに注意 */
 		XX = XX - 1; A = 0;                /* 連分数展開 */
 		for(i = N; i>0;i--)  
			A = (double)i * XX / (2 + (double)i * XX / (2 *(double) i + 1.0 + A));
 		return L + XX / (1.0 + A);
	}
(2)
	double mySqrt(double X){
   		double DX, DN, DM, Z; int i, N, NN;
   		double XX = fabs(X); if( XX < 1E-15) return 0.0;
   		bool less1 = false;
   		if( XX < 1){ XX = 1 / XX; less1 = true;};       /* 1より小さい時，逆数　　*/
   		for(N = 1, NN = 1; NN < XX; N++, NN = N * N);   /* 平方がXXを超えない整数 */
   		if(fabs((double)NN - XX) < 1E-15) return (double) N;
   		N--; NN = N * N; DX = XX - (double)NN;
   		if(fabs(DX) < 1E-15)              return (double) N;
   		DN = 2 * N; DM = DN / DX;                       /* 連分数の計算 */
   		for(Z = 0, i = 0; i < 20; i++) Z = 1 / (DM + 1 / (DN + Z));
   		Z = N + Z;                        /* 元データが1より小さい時，逆数をとる */
   		if(less1) return 1 / Z; else return Z;
	}

2.7	πを求める
(1)
	double getPi1(){
		int K ; double P = 0, bP, T, D;
		K = 1; T = 16.0 / 5.0; D = -5.0 * 5.0;
		do{ bP=P; P += T / K; T /= D; K += 2; } while(P != bP);
		K =1 ; T = 4.0/239.0; D = -239.0*239.0;
		do{ bP=P; P -= T / K; T /= D; K += 2; } while(P != bP);
		return P;
	}
(2)
	double getPi2(){
		int NN=30; double P = (double)(NN * 2 - 1);
		for(int i=NN; i>=1; i--) P = 2.0*(double)i -1.0+(double)(i * i) / P;
		return 4/P;
	}

2.8	ネイピア数(自然対数の底) eを求める
(1)
	double getE1(){
		double E, A, bE;int N; E=0; A=1; N=1;
		do{ bE = E; E +=A; A /= N; N++;} while(E != bE);
		return E;
	}
(2)
	double getE2(){
		double E = 1; int NN=20;
		for(int i = NN; i >= 2; i -= 2) E = 1.0 + 1.0 /((double)I + 1.0/(1.0 + 1.0/E));
		return 2 + 1 / E;
	}

2.9	配列の加算
(1)
	float sum1(int N,float A[]){
		float S=0;
		for(int i=0; i < N; i++)   S += A[i];
		return S;
	}
(2)
	float sum3(int N,float A[]){
		float S=0, R=0, T; 	/* 加算結果と積み残し分を0にする	*/
		for(int i=0; i <N; i++){
	  		R += A[i] ; 	/* 加算対象と積み残し分を加える	*/
	  		T  = S    ; 	/* それまでの和を保存	*/
	  		S += R    ; 	/* Rを加算したのに，実際には	*/
	  		T  = S - T; 	/* Tしか加算されていないので，	*/
	  		R -= T    ; 	/* 積み残し分(R - T)を計算する	*/
		}	
    		return S;
	}

2.10	ラグランジェ補間

	double LagInt(double X, int N, double Xtable[],double Ytable[]){
		double R = 0;
		for(int i = 0; i < N; i++){
	  		double T = Ytable[i];
	  		for(int j = 0;j < N; j++) if(i != j) T *=(X - Xtable[j]) / (Xtable[i] - Xtable[j]);
	  		R += T;
		}
		return R;
	}
2.11	スプライン補間

【呼出し形式】
	cmpSpline ( 始点直前境界条件, 始点直前固定値, 終点直後境界条件, 
	終点直後固定値, 補間X値 , テーブルの大きさ, テーブル)	
	【境界条件】		
	�@	True	始点の直前の値または終点の直後の値に固定値を使用
	�A	False	微分値一定。（始点／終点付近は線形補間となる）

	double splineInt(int startB,double startV, int endB, double endV,
 		double X,int N,double table[]){
  		int ID = (int) X, NM = N - 1;
  		double F1, F2, F3, F4, DF, DF1, DF2,A1, A2, A3, A4, H;
  		if(ID < 0) ID = 0; if(ID > NM) ID = NM;
  		F2 = table[ID];
  		if(ID == 0){
			F1 = startV; F3 = table[ID + 1]; F4 = table[ID + 2];
	        	if(!startB) F1 = 2 * F2 - F3;
  		}
  		else{
	  		F1 = table[ID - 1]; F3 = endV; F4 = endV;
	  		If     (ID == NM)      { 
				if(!endB){
					DF = F1 - F2; F3 = F2 - DF; F4 = F3 - DF;
				}
			}
	 		else if(ID == (NM - 1)){ F3 = table[ID + 1]; 
			if(!endB){ DF = F2 - F3; F4 = F3 - DF; }
	  		else     { F3 = table[ID + 1]; F4 = table[ID + 2];}
 		}
 		DF1 = (F3 - F1) / 2; DF2 = (F4 - F2) / 2; H = X - (double)ID;
  		A1  = 2 * (F2 - F3) + (DF1 + DF2);
  		A2  = 3 * (F3 - F2) - (2 * DF1 + DF2);
 		A3  = DF1; A4 = F2;
  		return  ((A1 * H + A2) * H + A3) * H + A4;
	}