ADMB Documentation  11.6rc.3404
 All Classes Files Functions Variables Typedefs Friends Defines
df3fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 
00008 #include <df1b2fun.h>
00009 //#include "df3fun.h"
00010   df1b2variable * df3_one_variable::ind_var=0;
00011 
00015 df3_one_variable::df3_one_variable()
00016 {
00017   v[0] = 0;
00018   v[1] = 0;
00019   v[2] = 0;
00020   v[3] = 0;
00021 }
00025 df3_one_variable::df3_one_variable(const df3_one_variable& x)
00026 {
00027   v[0] = x.v[0];
00028   v[1] = x.v[1];
00029   v[2] = x.v[2];
00030   v[3] = x.v[3];
00031 }
00035 df3_one_vector::df3_one_vector(const df3_one_vector& m2)
00036 {
00037   index_min = m2.index_min;
00038   index_max = m2.index_max;
00039   shape = m2.shape;
00040   if (shape)
00041   {
00042     (shape->ncopies)++;
00043   }
00044   v = m2.v;
00045 }
00049 df3_one_vector::~df3_one_vector()
00050 {
00051   deallocate();
00052 }
00053 
00054  void df3_one_vector::deallocate(void)
00055  {
00056    if(shape)
00057    {
00058      if (shape->ncopies)
00059      {
00060        (shape->ncopies)--;
00061      }
00062      else
00063      {
00064        v = (df3_one_variable*) (shape->trueptr);
00065        delete [] v;
00066        v = NULL;
00067        delete shape;
00068        shape=0;
00069      }
00070    }
00071  }
00072 
00073  dvector value(const df3_one_vector& v)
00074  {
00075    int mmin=v.indexmin();
00076    int mmax=v.indexmax();
00077    dvector cv(mmin,mmax);
00078    for (int i=mmin;i<=mmax;i++)
00079    {
00080      cv(i)=value(v(i));
00081    }
00082    return cv;
00083  }
00084 
00085  dvector first_derivatives(const df3_one_vector& v)
00086  {
00087    int mmin=v.indexmin();
00088    int mmax=v.indexmax();
00089    dvector cv(mmin,mmax);
00090    for (int i=mmin;i<=mmax;i++)
00091    {
00092      cv(i)=*(v(i).get_udot());
00093    }
00094    return cv;
00095  }
00096 
00097  dvector second_derivatives(const df3_one_vector& v)
00098  {
00099    int mmin=v.indexmin();
00100    int mmax=v.indexmax();
00101    dvector cv(mmin,mmax);
00102    for (int i=mmin;i<=mmax;i++)
00103    {
00104      cv(i)=*(v(i).get_udot2());
00105    }
00106    return cv;
00107  }
00108 
00109  dvector third_derivatives(const df3_one_vector& v)
00110  {
00111    int mmin=v.indexmin();
00112    int mmax=v.indexmax();
00113    dvector cv(mmin,mmax);
00114    for (int i=mmin;i<=mmax;i++)
00115    {
00116      cv(i)=*(v(i).get_udot3());
00117    }
00118    return cv;
00119  }
00120 
00121   void df3_one_vector::initialize(void)
00122   {
00123     int mmin=indexmin();
00124     int mmax=indexmax();
00125     for (int i=mmin;i<=mmax;i++)
00126     {
00127       (*this)(i)=0.0;
00128     }
00129   }
00130 
00131   df3_one_vector::df3_one_vector(void)
00132   {
00133     allocate();
00134   }
00135 
00136   df3_one_vector::df3_one_vector(int min,int max)
00137   {
00138     allocate(min,max);
00139   }
00140 
00141   void df3_one_vector::allocate(int min,int max)
00142   {
00143     index_min=min;
00144     index_max=max;
00145     v=new df3_one_variable[max-min+1];
00146     if (v==0)
00147     {
00148       cerr << "error allocating memory in df3_one_vector" << endl;
00149       ad_exit(1);
00150     }
00151     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00152     {
00153       cerr << "Error trying to allocate memory for df3_one_vector"
00154            << endl;;
00155       ad_exit(1);
00156     }
00157     v-=min;
00158   }
00159 
00160   void df3_one_vector::allocate(void)
00161   {
00162     index_min=0;
00163     index_max=-1;
00164     v=0;
00165     shape=0;
00166   }
00167 
00168  dmatrix value(const df3_one_matrix& v)
00169  {
00170    int rmin=v.indexmin();
00171    int rmax=v.indexmax();
00172    dmatrix cm(rmin,rmax);
00173    for (int i=rmin;i<=rmax;i++)
00174    {
00175      int cmin=v(i).indexmin();
00176      int cmax=v(i).indexmax();
00177      cm(i).allocate(cmin,cmax);
00178      for (int j=cmin;j<=cmax;j++)
00179      {
00180        cm(i,j)=value(v(i,j));
00181      }
00182    }
00183    return cm;
00184  }
00185 
00186  dmatrix first_derivatives(const df3_one_matrix& v)
00187  {
00188    int rmin=v.indexmin();
00189    int rmax=v.indexmax();
00190    dmatrix cm(rmin,rmax);
00191    for (int i=rmin;i<=rmax;i++)
00192    {
00193      int cmin=v(i).indexmin();
00194      int cmax=v(i).indexmax();
00195      cm(i).allocate(cmin,cmax);
00196      cm(i)=first_derivatives(v(i));
00197    }
00198    return cm;
00199  }
00200 
00201  dmatrix second_derivatives(const df3_one_matrix& v)
00202  {
00203    int rmin=v.indexmin();
00204    int rmax=v.indexmax();
00205    dmatrix cm(rmin,rmax);
00206    for (int i=rmin;i<=rmax;i++)
00207    {
00208      int cmin=v(i).indexmin();
00209      int cmax=v(i).indexmax();
00210      cm(i).allocate(cmin,cmax);
00211      cm(i)=second_derivatives(v(i));
00212    }
00213    return cm;
00214  }
00215 
00216  dmatrix third_derivatives(const df3_one_matrix& v)
00217  {
00218    int rmin=v.indexmin();
00219    int rmax=v.indexmax();
00220    dmatrix cm(rmin,rmax);
00221    for (int i=rmin;i<=rmax;i++)
00222    {
00223      int cmin=v(i).indexmin();
00224      int cmax=v(i).indexmax();
00225      cm(i).allocate(cmin,cmax);
00226      cm(i)=third_derivatives(v(i));
00227    }
00228    return cm;
00229  }
00230 
00231 
00232  df3_one_matrix::df3_one_matrix(const df3_one_matrix& m2)
00233  {
00234    index_min=m2.index_min;
00235    index_max=m2.index_max;
00236    shape=m2.shape;
00237    if (shape)
00238    {
00239      (shape->ncopies)++;
00240    }
00241    v = m2.v;
00242  }
00243 
00244  df3_one_matrix::~df3_one_matrix()
00245  {
00246    deallocate();
00247  }
00248 
00249  void df3_one_matrix::deallocate(void)
00250  {
00251    if (shape)
00252    {
00253      if (shape->ncopies)
00254      {
00255        (shape->ncopies)--;
00256      }
00257      else
00258      {
00259        v = (df3_one_vector*) (shape->get_pointer());
00260        delete [] v;
00261        v=0;
00262        delete shape;
00263        shape=0;
00264      }
00265    }
00266  }
00267 
00268 
00269   void df3_one_matrix::initialize(void)
00270   {
00271     int mmin=indexmin();
00272     int mmax=indexmax();
00273     for (int i=mmin;i<=mmax;i++)
00274     {
00275       (*this)(i).initialize();
00276     }
00277   }
00278 
00279 
00280   df3_one_matrix::df3_one_matrix(int rmin,int rmax,int cmin,int cmax)
00281   {
00282     index_min=rmin;
00283     index_max=rmax;
00284     v=new df3_one_vector[rmax-rmin+1];
00285     if (v==0)
00286     {
00287       cerr << "error allocating memory in df3_one_matrix" << endl;
00288       ad_exit(1);
00289     }
00290     if ( (shape=new mat_shapex(v)) == NULL)
00291     {
00292       cerr << "Error trying to allocate memory for df3_one_vector"
00293            << endl;;
00294     }
00295     v-=rmin;
00296 
00297     for (int i=rmin;i<=rmax;i++)
00298     {
00299       v[i].allocate(cmin,cmax);
00300     }
00301   }
00302 
00303 /*
00304   df3_one_variable operator F(const df3_one_variable& x)
00305   {
00306     df3_one_variable z;
00307 
00308     *z.get_u() = ::F(*x.get_u());
00309 
00310     *z.get_udot() = ::D1F(*x.get_u())* *x.get_udot();
00311 
00312     *z.get_udot2() = ::D2F(*x.get_u())* square(*x.get_udot())
00313                    + ::D1F(*x.get_u())* *x.get_udot2();
00314 
00315     *z.get_udot3() = ::D3F(*x.get_u()) * cube(*x.get_udot())
00316                    + 3.0 * ::D2F(*x.get_u()) * *x.get_udot() * *x.get_udot2()
00317                    + ::D1F(*x.get_u()) * *x.get_udot3();
00318     return z;
00319   }
00320 
00321 */
00322 
00326 df3_one_variable& df3_one_variable::operator-=(double d)
00327 {
00328   *get_u() -= d;
00329   return *this;
00330 }
00334 df3_one_variable& df3_one_variable::operator-=(const df3_one_variable& _v)
00335 {
00336   *get_u() -= *_v.get_u();
00337   *get_udot() -= *_v.get_udot();
00338   *get_udot2() -= *_v.get_udot2();
00339   *get_udot3() -= *_v.get_udot3();
00340   return *this;
00341 }
00342 
00343 df3_one_variable operator-(const df3_one_variable& v)
00344 {
00345   df3_one_variable z;
00346   *z.get_u() = -(*v.get_u());
00347 
00348   *z.get_udot() = -(*v.get_udot());
00349   *z.get_udot2() = -(*v.get_udot2());
00350   *z.get_udot3() = -(*v.get_udot3());
00351   return z;
00352 }
00356 df3_one_variable& df3_one_variable::operator+=(const df3_one_variable& _v)
00357 {
00358   *get_u() += *_v.get_u();
00359   *get_udot() += *_v.get_udot();
00360   *get_udot2() += *_v.get_udot2();
00361   *get_udot3() += *_v.get_udot3();
00362   return *this;
00363 }
00364 
00365   df3_one_variable sqrt(const df3_one_variable& x)
00366   {
00367     df3_one_variable z;
00368     double u=::sqrt(*x.get_u());
00369     double xinv=1.0/(*x.get_u());
00370     double zp=0.5/u;
00371     double zp2=-0.5*zp*xinv;
00372     double zp3=-1.5*zp2*xinv;
00373 
00374     *z.get_u() = u;
00375 
00376     *z.get_udot() = zp* *x.get_udot();
00377 
00378     *z.get_udot2() = zp2 * square(*x.get_udot())
00379                    + zp * *x.get_udot2();
00380 
00381     *z.get_udot3() = zp3 * cube(*x.get_udot())
00382                    + 3.0 * zp2 * *x.get_udot() * *x.get_udot2()
00383                    + zp * *x.get_udot3();
00384 
00385     return z;
00386   }
00387 
00388 
00389 
00390   df3_one_variable exp(const df3_one_variable& x)
00391   {
00392     df3_one_variable z;
00393 
00394     *z.get_u() = ::exp(*x.get_u());
00395 
00396     *z.get_udot() = ::exp(*x.get_u())* *x.get_udot();
00397 
00398     *z.get_udot2() = ::exp(*x.get_u())* square(*x.get_udot())
00399                    + ::exp(*x.get_u())* *x.get_udot2();
00400 
00401     *z.get_udot3() = ::exp(*x.get_u()) * cube(*x.get_udot())
00402                    + 3.0 * ::exp(*x.get_u()) * *x.get_udot() * *x.get_udot2()
00403                    + ::exp(*x.get_u()) * *x.get_udot3();
00404     return z;
00405   }
00406 
00407   df3_one_variable log(const df3_one_variable& x)
00408   {
00409     df3_one_variable z;
00410 
00411     double xp=1.0/(*x.get_u());
00412     double xp2=-square(xp);
00413     double xp3=-2.0*xp*xp2;
00414 
00415     *z.get_u() = ::log(*x.get_u());
00416 
00417     *z.get_udot() = xp * *x.get_udot();
00418 
00419     *z.get_udot2() = xp2* square(*x.get_udot())
00420                    + xp * *x.get_udot2();
00421 
00422     *z.get_udot3() = xp3 * cube(*x.get_udot())
00423                    + 3.0 * xp2 * *x.get_udot() * *x.get_udot2()
00424                    + xp * *x.get_udot3();
00425     return z;
00426   }
00427 
00428   df3_one_variable inv(const df3_one_variable& x)
00429   {
00430     df3_one_variable z;
00431     double xinv=1.0/(*x.get_u());
00432     double zp=-xinv*xinv;
00433     double zp2=-2.0*zp*xinv;
00434     double zp3=-3.0*zp2*xinv;
00435 
00436     *z.get_u() = xinv;
00437     *z.get_udot() = zp * *x.get_udot();
00438 
00439     *z.get_udot2() = zp2 * square(*x.get_udot())
00440                    + zp * *x.get_udot2();
00441 
00442     *z.get_udot3() = zp3 * cube(*x.get_udot())
00443                    + 3.0 * zp2 * *x.get_udot() * *x.get_udot2()
00444                    + zp * *x.get_udot3();
00445 
00446     return z;
00447   }
00448 
00449   df3_one_variable& df3_one_variable::operator = (const df3_one_variable& x)
00450   {
00451     *get_u() = *x.get_u();
00452     *get_udot() = *x.get_udot();
00453     *get_udot2() = *x.get_udot2();
00454     *get_udot3() = *x.get_udot3();
00455     return *this;
00456   }
00457 
00458   df3_one_variable& df3_one_variable::operator = (double x)
00459   {
00460     *get_u() = x;
00461     *get_udot() = 0.0;
00462     *get_udot2() = 0.0;
00463     *get_udot3() = 0.0;
00464     return *this;
00465   }
00466 
00467 
00468   df3_one_variable operator * (const df3_one_variable& x,
00469     const df3_one_variable& y)
00470   {
00471     df3_one_variable z;
00472     *z.get_u() = *x.get_u() *  *y.get_u();
00473     *z.get_udot() = *x.get_u() * *y.get_udot()
00474                  + *x.get_udot() * *y.get_u();
00475 
00476     *z.get_udot2() = *x.get_udot2() * *y.get_u()
00477                    + 2.0 * *x.get_udot() * *y.get_udot()
00478                    +  *x.get_u() * *y.get_udot2();
00479 
00480     *z.get_udot3() = *x.get_udot3() * *y.get_u()
00481                    + 3.0 * *x.get_udot2() * *y.get_udot()
00482                    + 3.0 * *x.get_udot() * *y.get_udot2()
00483                    +  *x.get_u() * *y.get_udot3();
00484     return z;
00485   }
00486 
00487   df3_one_variable operator * (double x,
00488     const df3_one_variable& y)
00489   {
00490     df3_one_variable z;
00491     *z.get_u() = x *  *y.get_u();
00492     *z.get_udot() = x * *y.get_udot();
00493 
00494     *z.get_udot2() = x * *y.get_udot2();
00495 
00496     *z.get_udot3() = x * *y.get_udot3();
00497     return z;
00498   }
00499 
00500   df3_one_variable operator * (const df3_one_variable& x,
00501     double y)
00502   {
00503     df3_one_variable z;
00504     *z.get_u() = *x.get_u() *  y;
00505     *z.get_udot() =  *x.get_udot() * y;
00506 
00507     *z.get_udot2() = *x.get_udot2() * y;
00508 
00509     *z.get_udot3() = *x.get_udot3() * y;
00510     return z;
00511   }
00512 
00513 
00514   df3_one_variable operator / (const df3_one_variable& x,
00515     const df3_one_variable& y)
00516   {
00517     df3_one_variable u=inv(y);
00518     return x*u;
00519   }
00520 
00521   df3_one_variable operator / (double x,
00522     const df3_one_variable& y)
00523   {
00524     df3_one_variable u=inv(y);
00525     return x*u;
00526   }
00527 
00528   df3_one_variable operator + (double x,const df3_one_variable& y)
00529   {
00530     df3_one_variable z;
00531     *z.get_u() =  x + *y.get_u();
00532     *z.get_udot() =  *y.get_udot();
00533     *z.get_udot2() = *y.get_udot2();
00534     *z.get_udot3() = *y.get_udot3();
00535     return z;
00536   }
00537 
00538  //
00539  //   df3_one_variable operator / (const df3_one_variable& x,
00540  //     const df3_one_variable& y)
00541  //   {
00542  //     df3_one_variable z;
00543  //     double yinv =  1.0 / (*y.get_u());
00544  //     double yinv2 = yinv * yinv;
00545  //     double yinv3 = yinv * yinv2;
00546  //     doubl yd = *y.get_udot();
00547  //
00548  //     //*z.get_u() = *x.get_u() /  *y.get_u();
00549  //     *z.get_u() = *x.get_u() * yinv;
00550  //
00551  //     *z.get_udot() =  - (*x.get_u()) * yinv2 * yd
00552  //                   + *x.get_udot() * yinv;
00553  //
00554  //     *z.get_udot2() = *x.get_udot2() * yinv
00555  //                    - 2.0 * *x.get_udot() * yd * yinv2
00556  //                    + 2.0 * *x.get_u() * yinv3  * yd *yd
00557  //                    -  *x.get_u() * yinv2 * y.get_udot2();
00558  //
00559  //     *z.get_udot3() = *x.get_udot3() * yinv
00560  //                    + 3.0 * *x.get_udot2() * *y.get_udot()
00561  //                    + 3.0 * *x.get_udot() * *y.get_udot2()
00562  //                    +  *x.get_u() * *y.get_udot3();
00563  //   }
00564  //
00565 
00566 
00567   df3_one_variable operator + (const df3_one_variable& y,
00568     const double x)
00569   {
00570     df3_one_variable z;
00571     *z.get_u() =   x + *y.get_u();
00572     *z.get_udot() =  *y.get_udot();
00573     *z.get_udot2() = *y.get_udot2();
00574     *z.get_udot3() = *y.get_udot3();
00575     return z;
00576   }
00577 
00578   df3_one_variable operator + (const df3_one_variable& x,
00579     const df3_one_variable& y)
00580   {
00581     df3_one_variable z;
00582     *z.get_u() = *x.get_u() + *y.get_u();
00583     *z.get_udot() = *x.get_udot() + *y.get_udot();
00584     *z.get_udot2() = *x.get_udot2() + *y.get_udot2();
00585     *z.get_udot3() = *x.get_udot3() + *y.get_udot3();
00586     return z;
00587   }
00588 
00589   df3_one_variable operator - (const df3_one_variable& x,
00590     const df3_one_variable& y)
00591   {
00592     df3_one_variable z;
00593     *z.get_u() = *x.get_u() - *y.get_u();
00594     *z.get_udot() = *x.get_udot() - *y.get_udot();
00595     *z.get_udot2() = *x.get_udot2() - *y.get_udot2();
00596     *z.get_udot3() = *x.get_udot3() - *y.get_udot3();
00597     return z;
00598   }
00599 
00600   df3_one_variable operator - (const df3_one_variable& x,
00601     double y)
00602   {
00603     df3_one_variable z;
00604     *z.get_u() = *x.get_u() - y;
00605     *z.get_udot() = *x.get_udot();
00606     *z.get_udot2() = *x.get_udot2();
00607     *z.get_udot3() = *x.get_udot3();
00608     return z;
00609   }
00610 
00611   init_df3_one_variable::init_df3_one_variable(const df1b2variable& _v)
00612   {
00613     ADUNCONST(df1b2variable,v)
00614    /*
00615     if (num_ind_var>0)
00616     {
00617       cerr << "can only have 1 independent_variables in df3_one_variable"
00618        " function" << endl;
00619       ad_exit(1);
00620    }
00621    */
00622    ind_var=&v;
00623     *get_u() =  *v.get_u();
00624     *get_udot() = 1.0;
00625     *get_udot2() = 0.0;
00626     *get_udot3() = 0.0;
00627   }
00628 
00629   init_df3_one_variable::init_df3_one_variable(double v)
00630   {
00631     *get_u() =  v;
00632     *get_udot() = 1.0;
00633     *get_udot2() = 0.0;
00634     *get_udot3() = 0.0;
00635   }
00636 
00637 df3_one_matrix choleski_decomp(const df3_one_matrix& MM)
00638 {
00639   // kludge to deal with constantness
00640   df3_one_matrix & M= (df3_one_matrix &) MM;
00641   int rmin=M.indexmin();
00642   int cmin=M(rmin).indexmin();
00643   int rmax=M.indexmax();
00644   int cmax=M(rmin).indexmax();
00645   if (rmin !=1 || cmin !=1)
00646   {
00647     cerr << "minimum row and column inidices must equal 1 in "
00648       "df1b2matrix choleski_decomp(const df3_one_atrix& MM)"
00649          << endl;
00650     ad_exit(1);
00651   }
00652   if (rmax !=cmax)
00653   {
00654     cerr << "Error in df1b2matrix choleski_decomp(const df3_one_matrix& MM)"
00655       " Matrix not square" << endl;
00656     ad_exit(1);
00657   }
00658 
00659   int n=rmax-rmin+1;
00660   df3_one_matrix L(1,n,1,n);
00661 #ifndef SAFE_INITIALIZE
00662     L.initialize();
00663 #endif
00664 
00665   int i,j,k;
00666   df3_one_variable tmp;
00667 
00668     if (value(M(1,1))<=0)
00669     {
00670       cerr << "Error matrix not positive definite in choleski_decomp"
00671         <<endl;
00672       ad_exit(1);
00673     }
00674 
00675   L(1,1)=sqrt(M(1,1));
00676   for (i=2;i<=n;i++)
00677   {
00678     L(i,1)=M(i,1)/L(1,1);
00679   }
00680 
00681   for (i=2;i<=n;i++)
00682   {
00683     for (j=2;j<=i-1;j++)
00684     {
00685       tmp=M(i,j);
00686       for (k=1;k<=j-1;k++)
00687       {
00688         tmp-=L(i,k)*L(j,k);
00689       }
00690       L(i,j)=tmp/L(j,j);
00691     }
00692     tmp=M(i,i);
00693     for (k=1;k<=i-1;k++)
00694     {
00695       tmp-=L(i,k)*L(i,k);
00696     }
00697 
00698     if (value(tmp)<=0)
00699     {
00700       cerr << "Error matrix not positive definite in choleski_decomp"
00701         <<endl;
00702       ad_exit(1);
00703     }
00704 
00705     L(i,i)=sqrt(tmp);
00706   }
00707 
00708   return L;
00709 }
00710 
00711 df1b2matrix& df1b2matrix::operator = (const df3_one_matrix& M)
00712 {
00713   int rmin=M.indexmin();
00714   int rmax=M.indexmax();
00715   if (rmin != indexmin() || rmax != indexmax())
00716   {
00717     cerr << "unequal shape in "
00718      "df1b2matrix& df1b2matrix::operator = (const df3_one_matrix& M)"
00719       << endl;
00720     ad_exit(1);
00721   }
00722 
00723   for (int i=rmin;i<=rmax;i++)
00724   {
00725     (*this)(i)=M(i);
00726   }
00727   return *this;
00728 }
00729 
00730 df1b2vector& df1b2vector::operator = (const df3_one_vector& M)
00731 {
00732   int rmin=M.indexmin();
00733   int rmax=M.indexmax();
00734   if (rmin != indexmin() || rmax != indexmax())
00735   {
00736     cerr << "unequal shape in "
00737      "df1b2vector& df1b2vector::operator = (const df3_one_vector& M)"
00738       << endl;
00739     ad_exit(1);
00740   }
00741 
00742   for (int i=rmin;i<=rmax;i++)
00743   {
00744     (*this)(i)=M(i);
00745   }
00746   return *this;
00747 }
00748 
00749 df1b2variable& df1b2variable::operator = (const df3_one_variable& v)
00750 {
00751   df1b2variable * px=df3_one_variable::ind_var;
00752   //df3_one_variable::ind_var=0;
00753   //df1b2variable * px=0;
00754   double  df= *v.get_udot();
00755   double * xd=px->get_u_dot();
00756   double * zd=get_u_dot();
00757   *get_u()=*v.get_u();
00758   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00759   {
00760     *zd++ = df * *xd++;
00761   }
00762 
00763  /*
00764   cout << *v.get_u()  << " ";
00765   cout << *v.get_udot()  << " ";
00766   cout << *v.get_udot2()  << " ";
00767   cout << *v.get_udot3()  << endl;
00768   */
00769   f1b2gradlist->write_pass1(px,this,*v.get_udot(),*v.get_udot2(),
00770 
00771     *v.get_udot3());
00772   return *this;
00773 }
00774 
00775 
00776 df1b2variable cumd_norm(const df1b2variable& _x)
00777 {
00778   df1b2variable z;
00779   init_df3_one_variable x(_x);
00780 
00781   const double b1=0.319381530;
00782   const double b2=-0.356563782;
00783   const double b3=1.781477937;
00784   const double b4=-1.821255978;
00785   const double b5=1.330274429;
00786   const double p=.2316419;
00787 
00788   if (value(x)>=0)
00789   {
00790     df3_one_variable u1=p*x;
00791     df3_one_variable u2=1.+u1;
00792     df3_one_variable u=1./u2;
00793     df3_one_variable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00794     df3_one_variable tmp1=-0.3989422804*exp(-.5*x*x);
00795     z=1.0+tmp1*y;
00796   }
00797   else
00798   {
00799     df3_one_variable w=-x;
00800     df3_one_variable u=1./(1+p*w);
00801     df3_one_variable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00802     df3_one_variable tmp1=0.3989422804*exp(-.5*x*x);
00803     z=tmp1*y;
00804   }
00805   return z;
00806 }
00807 
00808 df1b2variable bounded_cumd_norm(const df1b2variable& _x, double beta)
00809 {
00810   df1b2variable z;
00811   init_df3_one_variable x(_x);
00812 
00813   const double b1=0.319381530;
00814   const double b2=-0.356563782;
00815   const double b3=1.781477937;
00816   const double b4=-1.821255978;
00817   const double b5=1.330274429;
00818   const double p=.2316419;
00819 
00820   if (value(x)>=0)
00821   {
00822     df3_one_variable u1=p*x;
00823     df3_one_variable u2=1.+u1;
00824     df3_one_variable u=1./u2;
00825     df3_one_variable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00826     df3_one_variable tmp1=-0.3989422804*exp(-.5*x*x);
00827     z=1.0+tmp1*y;
00828   }
00829   else
00830   {
00831     df3_one_variable w=-x;
00832     df3_one_variable u=1./(1+p*w);
00833     df3_one_variable y=  ((((b5*u+b4)*u+b3)*u+b2)*u+b1)*u;
00834     df3_one_variable tmp1=0.3989422804*exp(-.5*x*x);
00835     z=tmp1*y;
00836   }
00837   z=beta*(z-0.5)+0.5;
00838   return z;
00839 }
00840 
00841