複素平面の‘上’で５次の多項式を描いたらどんなだろう

まず

結果
こんな感じ

で，

動機
タイトルのままの動機．ベジエ曲線と点との距離で５次方程式なんてのを見ていて思ったこと．

方針
とりあえずGSLで求めて（サンプリングして）gnuplotで描いてみるとかやってたんだけどピンとこなかったので，Quartz ComposerでOpenGLでポリゴン描かせてみることにした．（単なる複素数の積とか和の計算なので別にGSLを用いる必要はない．前回利用したのでその延長，また参考のサンプルを利用するためだけ）
gnuplotで描く時は求まった値の実部と虚部でベクトル場表示とかしていたのだがとりあえず今回は絶対値をy軸（高さ）にしてみる．そのままだとあまりに大きくなるので対数もとっている．
（x軸（横軸）が与えるxの実部，z軸（奥行き）が与えるxの虚部）

参考 ありがとうございました．
OpenGL プログラミング -- 3D 関数グラフィックス sample4
GLUTでOpenGL 3Dプログラミング | 法線ベクトルを計算する方法
（GSLのサンプル）

やったこと
5*x^5+4*x^4+3*x^3+2*x^2+1*x+0 の面を描く（たぶん）．

...
#include <math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
...
@implementation AlgebraicEquationPlugIn (Add)
GLfloat v[2000][2000][3];
GLfloat u[2000][2000][3];
gsl_complex _gsl_poly_complex_eval(double coef[], const int len, const gsl_complex x)
{
 gsl_complex ret;
 int i;
 
 /* Horner's Method */
 GSL_SET_COMPLEX(&ret, coef[len - 1], 0.0);
 
 /* ret = ret * x + coef[i] */
 for(i = len - 2; i >= 0; i--)
  ret = gsl_complex_add_real(gsl_complex_mul(ret, x), coef[i]);
 
 return ret;
}
int CalculateNormal(const float p1[3],
     const float p2[3],
     const float p3[3],     
     float n[3])
{
 float v1[3];
 float v2[3];
 float cross[3];
 float length;
 int i;
 
 /* v1 = p1 - p2を求める */
 for (i = 0; i < 3; i++) {
  v1[i] = p1[i] - p2[i];
 }
 
 /* v2 = p3 - p2を求める */
 for (i = 0; i < 3; i++) {
  v2[i] = p3[i] - p2[i];
 }
 
 /* 外積v2×v1（= cross）を求める */
 for (i = 0; i < 3; i++) {
  cross[i] = v2[(i+1)%3] * v1[(i+2)%3] - v2[(i+2)%3] * v1[(i+1)%3];
 }
 
 /* 外積v2×v1の長さ|v2×v1|（= length）を求める */
 length = sqrtf(cross[0] * cross[0] + 
                cross[1] * cross[1] + 
                cross[2] * cross[2]);
 
 /* 長さ|v2×v1|が0のときは法線ベクトルは求められない */
 if (length == 0.0f) {
  return 0;
 }
 
 /* 外積v2×v1を長さ|v2×v1|で割って法線ベクトルnを求める */
 for (i = 0; i < 3; i++) {
  n[i] = cross[i] / length;
 }
 return 1; 
}
void compute(void)
{
 int n;
 int x, z;
 double coef[6];
 gsl_complex csol, eval;
 
 for(n= 0; n<6; n++)
  coef[n] = (double)n;
 
 for(x= 0; x < 2000; x++) {
  for(z= 0; z < 2000; z++) {
   double d_x= ((double)x/1000.0)-1.0;
   double d_z= ((double)z/1000.0)-1.0;
   GSL_SET_COMPLEX(&csol, d_x, d_z);
   eval = _gsl_poly_complex_eval(coef, 6, csol);
   double r= GSL_REAL(eval);
   double i= GSL_IMAG(eval);
   v[x][z][0] = (GLfloat)d_x; 
   v[x][z][1] = (GLfloat)(log(sqrt(r*r+i*i))); 
   v[x][z][2] = (GLfloat)d_z;
  }
 }
 for(x= 0; x < 1999; x++) {
  for(z= 0; z < 1999; z++) {
   CalculateNormal(v[x+1][z], v[x][z], v[x][z+1], u[x][z]);
  }
 } 
}
- (void)draw:(id<QCPlugInContext>)context
{
 CGLContextObj cgl_ctx = [context CGLContextObj]; 
 if(cgl_ctx == NULL) 
  return; 

 glDisable(GL_CULL_FACE);
 glEnable(GL_DEPTH_TEST);
 GLfloat mat_ambient[]    = { 0.2f, 0.2f, 0.2f, 1.0f };
 GLfloat mat_diffuse[]    = { 0.8f, 0.8f, 1.0f, 1.0f };
 GLfloat mat_specular[]   = { 0.8f, 0.8f, 1.0f, 1.0f };
 GLfloat mat_shininess[]  = { 50 };
 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT,   mat_ambient); 
 glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE,   mat_diffuse);
 glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR,  mat_specular);
 glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess);
 
 int i,j;
 for(i= 0; i<1999; i++){
  glBegin(GL_QUAD_STRIP);
  for(j= 0; j<1999; j++){
   glNormal3f(u[i  ][j][0], u[i  ][j][1], u[i  ][j][2]); 
   glVertex3f(v[i  ][j][0], v[i  ][j][1], v[i  ][j][2]);
   glNormal3f(u[i+1][j][0], u[i+1][j][1], u[i+1][j][2]); 
   glVertex3f(v[i+1][j][0], v[i+1][j][1], v[i+1][j][2]);
  }
  glEnd();
 }
 glDisable(GL_DEPTH_TEST);
 
 return;
} 
@end
@implementation AlgebraicEquationPlugIn (Execution)
...
- (BOOL) startExecution:(id<QCPlugInContext>)context
{
 compute(); 
 return YES;
}
...
- (BOOL) execute:(id<QCPlugInContext>)context atTime:(NSTimeInterval)time withArguments:(NSDictionary*)arguments
{
 CGLContextObj cgl_ctx = [context CGLContextObj]; 
 if(cgl_ctx == NULL) 
  return NO;
 
 GLint saveMode; 
 glGetIntegerv(GL_MATRIX_MODE, &saveMode);
 glMatrixMode(GL_MODELVIEW); 
 glPushMatrix(); 
 
 [self draw:context];
 
 glMatrixMode(GL_MODELVIEW); 
 glPopMatrix(); 
 glMatrixMode(saveMode); 
 
 GLenum error;
 if(error = glGetError()) 
  [context logMessage:@"OpenGL error %04X", error];

 return (error ? NO : YES);  
}
...
@end

まぁそんだけ．