fork download
copy 
  1. #include <stdio.h>
  2. #include<math.h>
  3.  
  4.  
  5. typedef struct { double val[3]; } vect;
  6.  
  7.  
  8. typedef double (*Func)(vect);
  9.  
  10. #define EPSILON 0.00001
  11. #define N 3                    //未知数の数
  12.  
  13. int printA(double a[N][N+1]);
  14. void gauss(double a[N][N+1], double x[N]);
  15. int gauss_jordan(int n, double a[N][N+1], double b[N]);
  16.  
  17.  
  18. int printV(vect X){ for(int i=0;i<3;i++) printf("%10.3f",X.val[i]);	  printf("\n"); }
  19.  
  20.  
  21. double f1(vect X){ 
  22. double x=X.val[0],y=X.val[1],z=X.val[2];
  23. return -y*y * z*z - y*y + 24*y*z - z*z -13; }
  24.  
  25. double f2(vect X){ 
  26. double x=X.val[0],y=X.val[1],z=X.val[2];
  27. return -x*x * z*z - x*x + 24*x*z - z*z -13; }
  28.  
  29. double f3(vect X){ 
  30. double x=X.val[0],y=X.val[1],z=X.val[2];
  31. return -x*x * y*y - x*x + 24*x*y - y*y -13; }
  32.  
  33.  
  34. vect d2v(double a,double b,double c) { vect x; x.val[0]=a; x.val[1]=b; x.val[2]=c; return x; }
  35. vect vect_add(vect x, vect y) { vect z; z.val[0]=x.val[0]+y.val[0];  z.val[1]=x.val[1]+y.val[1];  z.val[2]=x.val[2]+y.val[2]; return z; }
  36. vect vect_hiku(vect x, vect y) { vect z; z.val[0]=x.val[0]-y.val[0];  z.val[1]=x.val[1]-y.val[1];  z.val[2]=x.val[2]-y.val[2]; return z; }
  37.  
  38. double bibun_x(Func f, vect a) { return  ( f(vect_add(a,d2v(EPSILON,0,0))) - f(a) )/EPSILON;  }
  39. double bibun_y(Func f, vect a) { return  ( f(vect_add(a,d2v(0,EPSILON,0))) - f(a) )/EPSILON;  }
  40. double bibun_z(Func f, vect a) { return  ( f(vect_add(a,d2v(0,0,EPSILON))) - f(a) )/EPSILON;  }
  41.  
  42.  
  43.  
  44.  
  45. int main() {
  46. vect V=d2v(1,1,1);
  47. double a[3][4], X[3];
  48.  
  49.     for (int i = 0; i < 10; i++) {
  50.   a[0][0] = bibun_x(f1,V);  a[0][1] = bibun_y(f1,V);  a[0][2] = bibun_z(f1,V);
  51.   a[1][0] = bibun_x(f2,V);  a[1][1] = bibun_x(f2,V);  a[1][2] = bibun_x(f2,V);
  52.   a[2][0] = bibun_x(f3,V);  a[2][1] = bibun_x(f3,V);  a[2][2] = bibun_x(f3,V);
  53.   a[0][3] = V.val[0];   a[1][3] = V.val[1];   a[2][3] = V.val[2];
  54.   X[0] = V.val[0];   X[1] = V.val[1];   X[2] = V.val[2];
  55. //printV(V);printA(a);gauss(a,X);
  56.  
  57.  gauss_jordan(3, a, X);
  58.  
  59. V = vect_hiku(V, d2v(X[0],X[1],X[2]));
  60. printV(V);
  61. }
  62. return 0;
  63. }
  64.  
  65.  
  66.  
  67.  
  68. //http://s...content-available-to-author-only...e.jp/wiki/wiki.cgi?page=%CF%A2%CE%A9%B0%EC%BC%A1%CA%FD%C4%F8%BC%B0%A4%CE%B2%F2%CB%A1
  69.  
  70.  
  71. void gauss(double a[N][N+1], double x[N]){
  72.         int i,j,k,l,pivot;
  73.         double p,q,m,tmp;
  74.  
  75. 	/*縦方向のループ*/
  76.         for(i=0;i<N;i++) {
  77. 	  m=0;
  78. 	  pivot=i;
  79.  
  80. 	  /*i列中でi行目より下にある要素内で最大を選ぶ*/
  81. 	  for(l=i;l<N;l++) {
  82. 	    /*i列の中で一番値が大きい行を選ぶ*/
  83. 	    if(fabs(a[l][i])>m) {
  84. 	      m=fabs(a[l][i]);
  85. 	      pivot=l;
  86. 	    }
  87. 	  }
  88.  
  89. 	  /*pivotがiと違えば、行の入れ替え(これを部分ピボットと呼ぶ)*/
  90. 	  if(pivot!=i) {  
  91. 	    for(j=0;j<N+1;j++) {
  92. 	      tmp=a[i][j];        
  93. 	      a[i][j]=a[pivot][j];
  94. 	      a[pivot][j]=tmp;
  95. 	    }
  96. 	  }
  97.         }
  98.  
  99.         for(k=0;k<N;k++) {
  100. 	  p=a[k][k];              //対角要素を保存
  101. 	  /*対角要素は1になることがわかっているので直接代入*/
  102. 	  a[k][k]=1;      
  103.  
  104. 	  /*k行目でk+1列より右にあるすべての要素を対角成分で割る*/
  105. 	  for(j=k+1;j<N+1;j++) {
  106. 	    a[k][j]/=p;
  107. 	  }
  108.  
  109. 	  for(i=k+1;i<N;i++) {
  110. 	    q=a[i][k];
  111.  
  112. 	    for(j=k+1;j<N+1;j++) {
  113. 	      a[i][j]-=q*a[k][j];
  114. 	    }
  115. 	    /*0となることがわかっているので直接代入*/
  116. 	    a[i][k]=0;
  117. 	  }
  118.         }
  119.  
  120. 	/*解の計算*/
  121.         for(i=N-1;i>=0;i--) {
  122. 	  x[i]=a[i][N];
  123. 	  for(j=N-1;j>i;j--) {
  124. 	    x[i]-=a[i][j]*x[j];
  125. 	  }
  126.         }
  127.  
  128. }
  129.  
  130.  
  131. int printA(double a[N][N+1]){
  132.         int i,j;
  133.         for(i=0;i<N;i++) {
  134. 	  for(j=0;j<N+1;j++) {
  135. 	    printf("%10.3f",a[i][j]);
  136. 	  }
  137. 	  printf("\n");
  138.  
  139.         }
  140. }
  141.  
  142.  
  143.  
  144.  
  145.  
  146. //http://w...content-available-to-author-only...0.jp/yamamoto/lecture/2006/5E/Linear_eauations/gaussj_html/node2.html#SECTION00021000000000000000
  147. int gauss_jordan(int n, double a[N][N+1], double b[N]) {
  148. 	   int ipv, i, j;
  149. 	   double inv_pivot, temp;
  150. 	   double big;
  151. 	   int pivot_row, row[30];
  152.  
  153. 	   for(ipv=1 ; ipv <= n ; ipv++){
  154.  
  155. 	      /* ---- 最大値探索 ---------------------------- */
  156. 	      big=0.0;
  157. 	      for(i=ipv ; i<=n ; i++){
  158. 	         if(fabs(a[i][ipv]) > big){
  159. 	            big = fabs(a[i][ipv]);
  160. 	            pivot_row = i;
  161. 	         }
  162. 	      }
  163. 	      if(big == 0.0){return 0;}	
  164. 	      row[ipv] = pivot_row;
  165.  
  166. 	      /* ---- 行の入れ替え -------------------------- */
  167. 	      if(ipv != pivot_row){
  168. 	         for(i=1 ; i<=n ; i++){
  169. 	            temp = a[ipv][i];
  170. 	            a[ipv][i] = a[pivot_row][i];
  171. 	            a[pivot_row][i] = temp;
  172. 	         }
  173. 	         temp = b[ipv];
  174. 	         b[ipv] = b[pivot_row];
  175. 	         b[pivot_row] = temp;
  176. 	      }
  177.  
  178. 	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
  179. 	      inv_pivot = 1.0/a[ipv][ipv];
  180. 	      a[ipv][ipv]=1.0;
  181. 	      for(j=1 ; j <= n ; j++){
  182. 	         a[ipv][j] *= inv_pivot;
  183. 	      }
  184. 	      b[ipv] *= inv_pivot;
  185.  
  186. 	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
  187. 	      for(i=1 ; i<=n ; i++){
  188. 	         if(i != ipv){
  189. 	            temp = a[i][ipv];
  190. 		    a[i][ipv]=0.0;
  191. 	            for(j=1 ; j<=n ; j++){
  192. 	               a[i][j] -= temp*a[ipv][j];
  193. 	            }
  194. 	            b[i] -= temp*b[ipv];
  195. 	         }
  196. 	      }
  197.  
  198. 	   }
  199.  
  200. 	   /* ---- 列の入れ替え(逆行列) -------------------------- */
  201. 	   for(j=n ; j>=1 ; j--){
  202. 	      if(j != row[j]){
  203. 	         for(i=1 ; i<=n ; i++){
  204. 	            temp = a[i][j];
  205. 	            a[i][j]=a[i][row[j]];
  206. 	            a[i][row[j]]=temp;
  207. 	         }
  208. 	      }
  209. 	   }
  210.  
  211. 	   return 1;
  212. 	}
Success #stdin #stdout 0.01s 2728KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
     0.000407174320635843350368991838208.000-20356191858780983151579627520.000
     0.000-165832362277029673580525978670004532868693673623389370056704.0002036515345749074135515201414365184.000
     0.000-27500372378380015283970316401357856949188082409879241601268973812467169238863681131234242502327726138036731137239285760.000-165832362277029673580525978670004532868693673623389370056704.000
     0.000-756270480949566529653565048268677523979737945769492451136667646343362153014653443451010825663342632770075504633194321905362395921981628743105764726910808931806420654351806899266899384642383627978002581410883678252440819360981354744381440.000-27500372378380015283970316401357856949188082409879241601268973812467169238863681131234242502327726138036731137239285760.000
     0.000-756270480949566529653565048268677523979737945769492451136667646343362153014653443451010825663342632770075504633194321905362395921981628743105764726910808931806420654351806899266899384642383627978002581410883678252440819360981354744381440.000       nan
     0.000       nan       nan
     0.000       nan       nan
     0.000       nan       nan
     0.000       nan       nan
     0.000       nan       nan
https://ideone.com/p3PZY
language:
C++ (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!