// -*- mode:C++ ; compile-command: "g++ -I.. -g -c isom.cc " -*- 
  #include "giacPCH.h"
  
  /*
   *  Copyright (C) 2001,14 R. De Graeve, Institut Fourier, 38402 St Martin d'Heres
   *
   *  This program is free software; you can redistribute it and/or modify
   *  it under the terms of the GNU General Public License as published by
   *  the Free Software Foundation; either version 3 of the License, or
   *  (at your option) any later version.
   *
   *  This program is distributed in the hope that it will be useful,
   *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   *  GNU General Public License for more details.
   *
   *  You should have received a copy of the GNU General Public License
   *  along with this program. If not, see <http://www.gnu.org/licenses/>.
   */
  
  using namespace std;
  #include "isom.h"
  #include "gen.h"
  #include "vecteur.h"
  #include "derive.h"
  #include "subst.h"
  #include "usual.h"
  #include "symbolic.h"
  #include "sym2poly.h"
  #include "giacintl.h"
  
  #ifndef NO_NAMESPACE_GIAC
  namespace giac {
  #endif // ndef NO_NAMESPACE_GIAC
  
    vecteur isom(const vecteur & M,GIAC_CONTEXT){
      gen errcode=checkanglemode(contextptr);
      if (is_undef(errcode)) return vecteur(1,errcode);
      int n;
      n=int(M.size());
      vecteur I;
      // for (int i=0;i<n;i++){
      //vecteur li(n);
      //li[i]=1;
      //I.push_back(li);
      //}
      I=midn(n);
      gen b;
      if (recursive_normal(mtran(M)*M,contextptr)==I){
        b=mdet(M,contextptr);
      } 
      else {
        b=gen(0);
      }
      if (is_zero(b)) { 
        return(makevecteur(b));
      }
      if (M==-I) {
        //on a une symetrie centrale
        return(makevecteur(cst_pi,b));
      }
      if ((n==2) && (b==1)) {
        //M est une rotation plane
        gen t;
        if (is_strictly_positive(M[1][0],contextptr)) {
  	t=acos(M[0][0],contextptr);
        }
        else {
  	t=-acos(M[0][0],contextptr);
        }
        return(makevecteur(t,b));
      }
      if (n==2) {
        //M est une symetrie par rapport a la droite de vecteur normal v[0]
        vecteur v;
        mker(addvecteur(M,I),v,contextptr);
        return(makevecteur(v[0],b));
      }
      if (M.size()==3) {
        vecteur v;
        //v[0] est l'axe de la rotation si b=1
        //v[0] est la normale au plan de symetrie si b=-1 ou v=R3 tout gen
        if (!mker(subvecteur(M,multvecteur(b,I)),v,contextptr))
  	return vecteur(1,gendimerr(contextptr));
        vecteur nn(3);
        nn[0]=v[0][0]; nn[1]=v[0][1]; nn[2]=v[0][2];
        vecteur w;
        if  (is_zero(b+1)) {
  	//w est _VECTose de vect qui engendre le plan de symetrie 
  	if (!mker(subvecteur(M,I),w,contextptr))
  	  return vecteur(1,gendimerr(contextptr));
        }
        else {
  	w=v;
        }
        //w contient les vecteurs invariants par M
        if (w.size()==2){
  	//on a une symetrie plane
  	return(makevecteur(nn,b));
        }
        if (w.size()!=2){
  	//on a une rotation(b=1) ou le produit rotation*symetrie (b=-1)
  	gen t;
  	vecteur u(3);
  	//2*cos(theta)+b=trace(M)
  	t=rdiv(M[0][0]+M[1][1]+M[2][2]-b,2,contextptr);
  	//on cherche le signe de theta (t=cos(theta)), u est orth a nn 
  	if ((nn[0] !=0)||(nn[1] !=0)) {
  	  u[0]=-nn[1];
  	  u[1]=nn[0];
  	  u[2]=0;
  	}
  	else {
  	  u[0]=1;
  	  u[1]=0;
  	  u[2]=0;
  	}
  	vecteur A;
  	A.push_back(nn);
  	A.push_back(u);
  	//A est la matrice formee par nn u M*u
  	A.push_back(multmatvecteur(M,u));
  	gen s;
  	s=mdet(A,contextptr);
  	//s est du signe de sin(theta) 
  	if (is_positive(s,contextptr)) { 
  	  //si s est >=0 l'angle est acos(t)
  	  return(makevecteur(nn,acos(t,contextptr),b));
  	}
  	else {
  	  //si s<0 l'angle est -acos(t)
  	  return(makevecteur(nn,-acos(t,contextptr),b));
  	}
        }
      }
      return 0;
    }
    static gen symb_isom(const gen & args){
      return symbolic(at_isom,args);
    }
    gen _isom(const gen & args,GIAC_CONTEXT){
      if ( args.type==_STRNG && args.subtype==-1) return  args;
      if (!ckmatrix(args))
        return symb_isom(args);
      return isom(*args._VECTptr,contextptr);
    }
    static const char _isom_s []="isom";
    static define_unary_function_eval (__isom,&_isom,_isom_s);
    define_unary_function_ptr5( at_isom ,alias_at_isom,&__isom,0,true);
  
    static int mkisom_teste(gen& n,int b, int & d1){
      int d;
      //d est la valeur de teste c'est la dim de l'espace du mkisom (d=2 ou >=3 )
      //d>3 pour faire des isometries orthogonales par rapport a un hyperplan
      if (n.type==_VECT){
        // n est un vecteur et d1 est la dimension du vecteur n 
        //si n n'est pas un vecteur il le devient! (cf else ...donc d1>=1)
        vecteur e=*(n._VECTptr);
        d1=int(e.size());
        if (d1>=3) {
  	d=d1;
        }
        else{
  	if (d1==1){
  	  if (b==1){d=2;}else {d=3;}
  	}
  	else{
  	  if (n[0].type==_VECT){d=3;}else{d=2;}
  	}
        }
      }
      else{
        //ds ce cas n n'est pas un vecteur: il le devient! 
        n=makevecteur(n);
        d1=1;
        if (b==1) {d=2;}else {d=3;}  
      }
      return(d);
    }
  
    vecteur mkisom(const gen & n_orig,int b,GIAC_CONTEXT) {
      checkanglemode(contextptr); 
      //n=les elements caracteristiques de l'isometrie et b=+1(rot) ou b=-1
      int d;
      int d1;
      gen n(n_orig);
      d=mkisom_teste(n,b,d1);
      //d=2 pour les isometries de R2 et d=3 pour celles de R3;
      //d1=dimension de n (n est devenu un vecteur)
      if (d==2) {
        if (b==1){
  	gen theta;
  	theta=n[0];
  	vecteur M2;
  	vecteur li(2); 
  	li[0]=cos(theta,contextptr); 
  	li[1]=-sin(theta,contextptr);
  	M2.push_back(li);
  	li[0]=sin(theta,contextptr);
  	li[1]=cos(theta,contextptr);
  	M2.push_back(li);
  	return(M2);
        }
        else {
  	//if (b==-1){
  	gen a=n[0];
  	gen b=n[1];
  	vecteur M2;
  	vecteur li(2); 
  	li[0]=rdiv(b*b-a*a,a*a+b*b,contextptr);
  	li[1]=-rdiv(gen(2)*a*b,a*a+b*b,contextptr);
  	M2.push_back(li);
  	li[0]=-rdiv(gen(2)*a*b,a*a+b*b,contextptr);
  	li[1]=-rdiv(b*b-a*a,a*a+b*b,contextptr);
  	M2.push_back(li);
  	return(M2);
  	//}
        }
      }  
      if (d>=3) {
        vecteur S;
        vecteur R;
        if (d1==1){
  	//on a une symetrie point
  	vecteur I;
  	I=midn(d);
  	return negvecteur(I);
        }
        vecteur nn(d);
        //on a une symetrie plan ou une rotation ou le produit rotation symetrie
        if (d1>=3){
  	//on a une symetrie plan et nn=n est normal au plan
  	//gen norme2=dotvecteur(n,n);
  	nn=*(n._VECTptr);
        }
        else {
  	//on a une rot d'axe nn=n[0] ou rot*sym d'axe nn=n[0] et plan orth a nn
  	nn=*(n[0])._VECTptr;
        } 
        vecteur ntn(d);
        //ntn est la matrice d*d egale a : nn*transpose(nn)/(nn*nn)
        gen norme2=dotvecteur(nn,nn);
        //nnn est le vecteur nn divise par sa norme au carre
        vecteur nnn=divvecteur(nn, norme2);
        for (int i=0;i<d;i++){
  	ntn[i]=multvecteur(nn[i],nnn);
        }
        //ntn[1]=multvecteur(nn[1],nnn);
        //ntn[2]=multvecteur(nn[2],nnn);
        if (b==-1){
  	//S=I3-2*nn*transpose(nn)/(nn*nn)
  	S=*(recursive_normal(midn(d)-multvecteur(2,ntn),contextptr)._VECTptr);
        }
        if (d1>=3) {return(S);}  
        if (d1==2){
  	//on a une rotation si b=1
  	vecteur A;
  	//A a comme vecteurs colonnes : produitvectoriel(nnn,ei)
  	vecteur li(3);
  	li[0]=0;
  	li[1]=-nnn[2];
  	li[2]=nnn[1];
  	A.push_back(li);
  	li[0]=nnn[2];
  	li[1]=0;
  	li[2]=-nnn[0];
  	A.push_back(li);
  	li[0]=-nnn[1];
  	li[1]=nnn[0];
  	li[2]=0;
  	A.push_back(li);
  	// A*sqrt(norme2) pour avoir la matrice prod vect(n,ei) avec norme(n)=1;
  	A=multvecteur(sqrt(norme2,contextptr),A);
  	//cout<<A<<endl;
  	gen theta=n[1];
  	R=*(recursive_normal(multvecteur(cos(theta,contextptr),midn(3))+multvecteur(1-cos(theta,contextptr),ntn)+multvecteur(sin(theta,contextptr),A),contextptr)._VECTptr); 
        }
        if (b==1) { 
  	return(R);
        }
        else {
  	return *(recursive_normal(mmult(S,R),0)._VECTptr); 
        }
      }
      return 0;
    }
    static gen symb_mkisom(const gen & args){
      return symbolic(at_mkisom,args);
    }
    gen _mkisom(const gen & args,GIAC_CONTEXT){
      if ( args.type==_STRNG && args.subtype==-1) return  args;
      if (args.type!=_VECT)
        return symb_mkisom(args);
      int s=int(args._VECTptr->size());
      if (s!=2)
        return gendimerr();
      if (args._VECTptr->back().type==_INT_){
        gen n=args._VECTptr->front();
        int b=args._VECTptr->back().val;
        return mkisom(n,b,contextptr);
      }
      return symb_mkisom(args);
    }
    static const char _mkisom_s []="mkisom";
    static define_unary_function_eval (__mkisom,&_mkisom,_mkisom_s);
    define_unary_function_ptr5( at_mkisom ,alias_at_mkisom,&__mkisom,0,true);
  
  
  
  #ifndef NO_NAMESPACE_GIAC
  } // namespace giac
  #endif // ndef NO_NAMESPACE_GIAC