fork download 
  1. // 3元
  2. // -y＾2 * z＾2 - y＾2 + 24*y*z - z＾2 -13 = 0
  3. // -x＾2 * z＾2 - x＾2 + 24*x*z - z＾2 -13 = 0
  4. // -x＾2 * y＾2 - x＾2 + 24*x*y - y＾2 -13 = 0
  5.  
  6. #define _CRT_SECURE_NO_WARNINGS
  7. #include <stdio.h>
  8.  
  9. double f1(double x, double y, double z)
  10. {
  11. 	return -y*y * z*z - y*y + 24*y*z - z*z -13;
  12. }
  13. double f2(double x, double y, double z)
  14. {
  15. 	return -x*x * z*z - x*x + 24*x*z - z*z -13;
  16. }
  17. double f3(double x, double y, double z)
  18. {
  19. 	return -x*x * y*y - x*x + 24*x*y - y*y -13;
  20. }
  21. double J11(double x, double y, double z)	// df1(x,y,z)/dx
  22. {
  23. 	return 0;
  24. }
  25. double J12(double x, double y, double z)	// df1(x,y,z)/dy
  26. {
  27. 	return -2*z*z*y - 2*y + 24*z;
  28. }
  29. double J13(double x, double y, double z)	// df1(x,y,z)/dz
  30. {
  31. 	return -2*y*y*z + 24*y - 2*z;
  32. }
  33. double J21(double x, double y, double z)	// df2(x,y,z)/dx
  34. {
  35. 	return -2*z*z*x - 2*x + 24*z;
  36. }
  37. double J22(double x, double y, double z)	// df2(x,y,z)/dy
  38. {
  39. 	return 0;
  40. }
  41. double J23(double x, double y, double z)	// df2(x,y,z)/dz
  42. {
  43. 	return -2*x*x*z + 24*x - 2*z;
  44. }
  45. double J31(double x, double y, double z)	// df3(x,y,z)/dx
  46. {
  47. 	return -2*y*y*x - 2*x + 24*y;
  48. }
  49. double J32(double x, double y, double z)	// df3(x,y,z)/dy
  50. {
  51. 	return -2*x*x*y + 24*x - 2*y;
  52. }
  53. double J33(double x, double y, double z)	// df3(x,y,z)/dz
  54. {
  55. 	return 0;
  56. }
  57.  
  58. void ucalc(double x, double y, double z, double *ux, double *uy, double *uz)
  59. {
  60. 	double	det;
  61.  
  62. 	det	= J11(x,y,z) * J22(x,y,z) * J33(x,y,z)
  63. 		+ J21(x,y,z) * J32(x,y,z) * J13(x,y,z)
  64. 		+ J31(x,y,z) * J12(x,y,z) * J23(x,y,z)
  65. 		- J11(x,y,z) * J32(x,y,z) * J23(x,y,z)
  66. 		- J31(x,y,z) * J22(x,y,z) * J13(x,y,z)
  67. 		- J21(x,y,z) * J12(x,y,z) * J33(x,y,z);
  68. 	*ux = 1/det * (
  69. 		(J22(x,y,z)*J33(x,y,z)-J23(x,y,z)*J32(x,y,z))*f1(x,y,z) +
  70. 		(J13(x,y,z)*J32(x,y,z)-J12(x,y,z)*J33(x,y,z))*f2(x,y,z) +
  71. 		(J12(x,y,z)*J23(x,y,z)-J13(x,y,z)*J22(x,y,z))*f3(x,y,z));
  72. 	*uy = 1/det * (
  73. 		(J23(x,y,z)*J31(x,y,z)-J21(x,y,z)*J33(x,y,z))*f1(x,y,z) +
  74. 		(J11(x,y,z)*J33(x,y,z)-J13(x,y,z)*J31(x,y,z))*f2(x,y,z) +
  75. 		(J13(x,y,z)*J21(x,y,z)-J11(x,y,z)*J23(x,y,z))*f3(x,y,z));
  76. 	*uz = 1/det * (
  77. 		(J21(x,y,z)*J32(x,y,z)-J22(x,y,z)*J31(x,y,z))*f1(x,y,z) +
  78. 		(J12(x,y,z)*J31(x,y,z)-J11(x,y,z)*J32(x,y,z))*f2(x,y,z) +
  79. 		(J11(x,y,z)*J22(x,y,z)-J12(x,y,z)*J21(x,y,z))*f3(x,y,z));
  80. }
  81.  
  82. int main()
  83. {
  84. 	double	x, y, z, ux, uy, uz;
  85. 	int	i;
  86.  
  87. 	while (1) {
  88. 		printf("x y z を入力(0 0 0で終了)\n");
  89. 		scanf("%lf%lf%lf", &x, &y, &z);
  90. 		if (x == 0 && y == 0 && z == 0) break;
  91. 		for (i = 0; i < 10; i++) {
  92. 			ucalc(x, y, z, &ux, &uy, &uz);
  93. 			x -= ux;
  94. 			y -= uy;
  95. 			z -= uz;
  96. //			printf("%d %.20f %.20f %.20f\n", i, x, y, z);
  97. 		}
  98. 		printf("x=%.20f y=%.20f z=%.20f\n", x, y, z);
  99. 		printf("f1=%.20f\n", f1(x, y, z));
  100. 		printf("f2=%.20f\n", f2(x, y, z));
  101. 		printf("f3=%.20f\n", f3(x, y, z));
  102. 	}
  103. 	return 0;
  104. }
  105.  
Success #stdin #stdout 0.01s 1680KB
comments (?)
stdin 
1 1 1
-1 -1 -1
5 5 0
0 5 -5
0 0 0
stdout 
x y z を入力(0 0 0で終了)
x=0.77954804507915753664 y=0.77954804507915753664 z=0.77954804507915753664
f1=-0.00000000000000016567
f2=-0.00000000000000016567
f3=-0.00000000000000016567
x y z を入力(0 0 0で終了)
x=-0.77954804507915753664 y=-0.77954804507915753664 z=-0.77954804507915753664
f1=-0.00000000000000016567
f2=-0.00000000000000016567
f3=-0.00000000000000016567
x y z を入力(0 0 0で終了)
x=4.62518160134424061170 y=4.62518160134423883534 z=0.33207309836566406558
f1=-0.00000000000000100440
f2=-0.00000000000000508968
f3=-0.00000000000003970609
x y z を入力(0 0 0で終了)
x=nan y=nan z=nan
f1=nan
f2=nan
f3=nan
x y z を入力(0 0 0で終了)
https://ideone.com/L4Fk7
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!