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