gluUnProjectを自前実装する

template<
  typename real_tC,
  typename real_tA,
  typename real_tB>
real_tC* mult_matrix44(real_tC* C, const real_tA* A, const real_tB* B) {
  return mult_matrixIJK(C, A, B, 4, 4, 4);
}
template<
  typename real_tC,
  typename real_tA,
  typename real_tB>
  real_tC* mult_matrix_vec41(real_tC* C, const real_tA* A, const real_tB* B) {
  return mult_matrixIJK(C, A, B, 4, 4, 1);
}

///////////////////////////////////////
// gluUnProject
// https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/gluUnProject.xml
//! @brief gluUnProjectの自前実装版
//! @param [in] winX 画面上のX座標
//! @param [in] winY 画面上のY座標
//! @param [in] winZ 画面上のZ座標
//! @param [in] model モデルビュー行列
//! @param [in] proj プロジェクション行列
//! @param [in] ビューポート x,y,width,height
//! @param [out] objX 空間座標のX
//! @param [out] objY 空間座標のY
//! @param [out] objZ 空間座標のZ
//! @retval true 計算成功
//! @retval false 計算失敗
template<typename real_t>
bool myUnProject(
  real_t winX,
  real_t winY,
  real_t winZ,
  const real_t* model,
  const real_t* proj,
  const int* view,
  real_t* objX,
  real_t* objY,
  real_t* objZ
) {

  real_t pm[16];
  real_t invpm[16];

  mult_matrix44(pm, proj, model);//pm = P×M
  bool ret = inverse_matrix44(invpm, pm);// pmの逆行列

  if (ret == false)
    return false;

  real_t v[4] = {
    2 * (winX - view[0]) / view[2] - 1.0 ,
    2 * (winY - view[1]) / view[3] - 1.0 ,
    2 * winZ - 1.0,
    1
  };

  real_t vd[4];
  mult_matrix_vec41(vd, invpm, v); // vd = invpm × v

  *objX = vd[0] / vd[3];
  *objY = vd[1] / vd[3];
  *objZ = vd[2] / vd[3];

  return true;
}

int WIDTH = 300;
int HEIGHT = 400;
void disp(void) {

  glViewport(0, 0, WIDTH, HEIGHT);


  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
  glLoadIdentity();
  gluPerspective(45, 0.5, 0.1, 2.5);

  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
  glTranslated(-1, 0, -2);
  glRotated(1.5, 1.0, 1.0, 1.0);

  ////作成した変換行列を取得する
  GLdouble  modeld[16];
  glGetDoublev(GL_MODELVIEW_MATRIX, modeld);

  GLdouble  projd[16];
  glGetDoublev(GL_PROJECTION_MATRIX, projd);

  GLint view[4];
  glGetIntegerv(GL_VIEWPORT, view);

  // 三次元座標（出力）
  double
    objx,
    objy,
    objz;

  // 二次元ピクセル座標（入力）
  double 
    winx = (rand() % 1000) / 100.0,
    winy = (rand() % 1000) / 100.0,
    winz = (rand() % 1000) / 100.0;

  gluUnProject(
    winx, winy, winz,
    modeld, projd, view,
    &objx, &objy, &objz);

  printf("glu %lf %lf %lf\n", objx, objy, objz);

  myUnProject(
    winx, winy, winz,
    modeld, projd, view,
    &objx, &objy, &objz
  );

  printf("my  %lf %lf %lf\n", objx, objy, objz);

  puts("-----------");

  glMatrixMode(GL_PROJECTION);
  glPopMatrix();

  glMatrixMode(GL_MODELVIEW);
  glPopMatrix();

}