GNU GENERAL PUBLIC LICENSE
		       Version 2, June 1991

 Copyright (C) 1989, 1991 Free Software Foundation, Inc.
                       59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 Everyone is permitted to copy and distribute verbatim copies
 of this license document, but changing it is not allowed.

			    Preamble

  The licenses for most software are designed to take away your
freedom to share and change it.  By contrast, the GNU General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users.  This
General Public License applies to most of the Free Software
Foundation's software and to any other program whose authors commit to
using it.  (Some other Free Software Foundation software is covered by
the GNU Library General Public License instead.)  You can apply it to
your programs, too.

  When we speak of free software, we are referring to freedom, not
price.  Our General Public Licenses are designed to make sure that you
have the freedom to distribute copies of free software (and charge for
this service if you wish), that you receive source code or can get it
if you want it, that you can change the software or use pieces of it
in new free programs; and that you know you can do these things.

  To protect your rights, we need to make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.

  For example, if you distribute copies of such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have.  You must make sure that they, too, receive or can get the
source code.  And you must show them these terms so they know their
rights.

  We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.

  Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software.  If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.

  Finally, any free program is threatened constantly by software
patents.  We wish to avoid the danger that redistributors of a free
program will individually obtain patent licenses, in effect making the
program proprietary.  To prevent this, we have made it clear that any
patent must be licensed for everyone's free use or not licensed at all.

  The precise terms and conditions for copying, distribution and
modification follow.


		    GNU GENERAL PUBLIC LICENSE
   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION

  0. This License applies to any program or other work which contains
a notice placed by the copyright holder saying it may be distributed
under the terms of this General Public License.  The "Program", below,
refers to any such program or work, and a "work based on the Program"
means either the Program or any derivative work under copyright law:
that is to say, a work containing the Program or a portion of it,
either verbatim or with modifications and/or translated into another
language.  (Hereinafter, translation is included without limitation in
the term "modification".)  Each licensee is addressed as "you".

Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope.  The act of
running the Program is not restricted, and the output from the Program
is covered only if its contents constitute a work based on the
Program (independent of having been made by running the Program).
Whether that is true depends on what the Program does.

  1. You may copy and distribute verbatim copies of the Program's
source code as you receive it, in any medium, provided that you
conspicuously and appropriately publish on each copy an appropriate
copyright notice and disclaimer of warranty; keep intact all the
notices that refer to this License and to the absence of any warranty;
and give any other recipients of the Program a copy of this License
along with the Program.

You may charge a fee for the physical act of transferring a copy, and
you may at your option offer warranty protection in exchange for a fee.

  2. You may modify your copy or copies of the Program or any portion
of it, thus forming a work based on the Program, and copy and
distribute such modifications or work under the terms of Section 1
above, provided that you also meet all of these conditions:

    a) You must cause the modified files to carry prominent notices
    stating that you changed the files and the date of any change.

    b) You must cause any work that you distribute or publish, that in
    whole or in part contains or is derived from the Program or any
    part thereof, to be licensed as a whole at no charge to all third
    parties under the terms of this License.

    c) If the modified program normally reads commands interactively
    when run, you must cause it, when started running for such
    interactive use in the most ordinary way, to print or display an
    announcement including an appropriate copyright notice and a
    notice that there is no warranty (or else, saying that you provide
    a warranty) and that users may redistribute the program under
    these conditions, and telling the user how to view a copy of this
    License.  (Exception: if the Program itself is interactive but
    does not normally print such an announcement, your work based on
    the Program is not required to print an announcement.)


These requirements apply to the modified work as a whole.  If
identifiable sections of that work are not derived from the Program,
and can be reasonably considered independent and separate works in
themselves, then this License, and its terms, do not apply to those
sections when you distribute them as separate works.  But when you
distribute the same sections as part of a whole which is a work based
on the Program, the distribution of the whole must be on the terms of
this License, whose permissions for other licensees extend to the
entire whole, and thus to each and every part regardless of who wrote it.

Thus, it is not the intent of this section to claim rights or contest
your rights to work written entirely by you; rather, the intent is to
exercise the right to control the distribution of derivative or
collective works based on the Program.

In addition, mere aggregation of another work not based on the Program
with the Program (or with a work based on the Program) on a volume of
a storage or distribution medium does not bring the other work under
the scope of this License.

  3. You may copy and distribute the Program (or a work based on it,
under Section 2) in object code or executable form under the terms of
Sections 1 and 2 above provided that you also do one of the following:

    a) Accompany it with the complete corresponding machine-readable
    source code, which must be distributed under the terms of Sections
    1 and 2 above on a medium customarily used for software interchange; or,

    b) Accompany it with a written offer, valid for at least three
    years, to give any third party, for a charge no more than your
    cost of physically performing source distribution, a complete
    machine-readable copy of the corresponding source code, to be
    distributed under the terms of Sections 1 and 2 above on a medium
    customarily used for software interchange; or,

    c) Accompany it with the information you received as to the offer
    to distribute corresponding source code.  (This alternative is
    allowed only for noncommercial distribution and only if you
    received the program in object code or executable form with such
    an offer, in accord with Subsection b above.)

The source code for a work means the preferred form of the work for
making modifications to it.  For an executable work, complete source
code means all the source code for all modules it contains, plus any
associated interface definition files, plus the scripts used to
control compilation and installation of the executable.  However, as a
special exception, the source code distributed need not include
anything that is normally distributed (in either source or binary
form) with the major components (compiler, kernel, and so on) of the
operating system on which the executable runs, unless that component
itself accompanies the executable.

If distribution of executable or object code is made by offering
access to copy from a designated place, then offering equivalent
access to copy the source code from the same place counts as
distribution of the source code, even though third parties are not
compelled to copy the source along with the object code.


  4. You may not copy, modify, sublicense, or distribute the Program
except as expressly provided under this License.  Any attempt
otherwise to copy, modify, sublicense or distribute the Program is
void, and will automatically terminate your rights under this License.
However, parties who have received copies, or rights, from you under
this License will not have their licenses terminated so long as such
parties remain in full compliance.

  5. You are not required to accept this License, since you have not
signed it.  However, nothing else grants you permission to modify or
distribute the Program or its derivative works.  These actions are
prohibited by law if you do not accept this License.  Therefore, by
modifying or distributing the Program (or any work based on the
Program), you indicate your acceptance of this License to do so, and
all its terms and conditions for copying, distributing or modifying
the Program or works based on it.

  6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the
original licensor to copy, distribute or modify the Program subject to
these terms and conditions.  You may not impose any further
restrictions on the recipients' exercise of the rights granted herein.
You are not responsible for enforcing compliance by third parties to
this License.

  7. If, as a consequence of a court judgment or allegation of patent
infringement or for any other reason (not limited to patent issues),
conditions are imposed on you (whether by court order, agreement or
otherwise) that contradict the conditions of this License, they do not
excuse you from the conditions of this License.  If you cannot
distribute so as to satisfy simultaneously your obligations under this
License and any other pertinent obligations, then as a consequence you
may not distribute the Program at all.  For example, if a patent
license would not permit royalty-free redistribution of the Program by
all those who receive copies directly or indirectly through you, then
the only way you could satisfy both it and this License would be to
refrain entirely from distribution of the Program.

If any portion of this section is held invalid or unenforceable under
any particular circumstance, the balance of the section is intended to
apply and the section as a whole is intended to apply in other
circumstances.

It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any
such claims; this section has the sole purpose of protecting the
integrity of the free software distribution system, which is
implemented by public license practices.  Many people have made
generous contributions to the wide range of software distributed
through that system in reliance on consistent application of that
system; it is up to the author/donor to decide if he or she is willing
to distribute software through any other system and a licensee cannot
impose that choice.

This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.


  8. If the distribution and/or use of the Program is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Program under this License
may add an explicit geographical distribution limitation excluding
those countries, so that distribution is permitted only in or among
countries not thus excluded.  In such case, this License incorporates
the limitation as if written in the body of this License.

  9. The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time.  Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.

Each version is given a distinguishing version number.  If the Program
specifies a version number of this License which applies to it and "any
later version", you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation.  If the Program does not specify a version number of
this License, you may choose any version ever published by the Free Software
Foundation.

  10. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission.  For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this.  Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.

			    NO WARRANTY

  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.

  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.

		     END OF TERMS AND CONDITIONS


	    How to Apply These Terms to Your New Programs

  If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
free software which everyone can redistribute and change under these terms.

  To do so, attach the following notices to the program.  It is safest
to attach them to the start of each source file to most effectively
convey the exclusion of warranty; and each file should have at least
the "copyright" line and a pointer to where the full notice is found.

    <one line to give the program's name and a brief idea of what it does.>
    Copyright (C) <year>  <name of author>

    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 2 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, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


Also add information on how to contact you by electronic and paper mail.

If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:

    Gnomovision version 69, Copyright (C) year name of author
    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
    This is free software, and you are welcome to redistribute it
    under certain conditions; type `show c' for details.

The hypothetical commands `show w' and `show c' should show the appropriate
parts of the General Public License.  Of course, the commands you use may
be called something other than `show w' and `show c'; they could even be
mouse-clicks or menu items--whatever suits your program.

You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary.  Here is a sample; alter the names:

  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
  `Gnomovision' (which makes passes at compilers) written by James Hacker.

  <signature of Ty Coon>, 1 April 1989
  Ty Coon, President of Vice

This General Public License does not permit incorporating your program into
proprietary programs.  If your program is a subroutine library, you may
consider it more useful to permit linking proprietary applications with the
library.  If this is what you want to do, use the GNU Library General
Public License instead of this License.

 
	./build_sources.sh

 
	./configure --prefix=/usr/local

 
	make

 
	make test

 
	make install

 

// standard input output library
# include <iostream>

// f2cad set up for adouble
# define F2CAD_USE_ADOLC

// prototype for dscal
# include <f2cad/dscal.hpp>

int main(void)
{	
	// turn on recording of adouble operations
	int tag  = 1;
	int keep = 1;
	trace_on(tag, keep);

	// independent variables
	int n     = 1;   // number of independent variables
	adouble  a[1];   // vector of independent variables
	a[0]   <<= 5.;   // first independent variables

	// create data structure expected by dscal
	integer  m = 2;
	adouble  f[2];
	f[0]     = 1.;
	f[1]     = 2.;
	integer incf = 1;

	// set f = a * f 
	f2cad::dscal_(&m, a, f, &incf );

	// dependent variables
	double fi;
	f[0] >>= fi;     // first dependent variable value
	f[1] >>= fi;     // second dependent variable value

	// turn off recording of adouble operations
	trace_off();

	// set differential for a 
	int d    = 1;    // order of the differential  
	double *daMemory = new double [2 * n];
	double **da      = new double * [n];
	int i, j;
	for(j = 0; j < n; j++)
		da[j] = daMemory + 2 * j;
	da[0][0] = .5;   // value of a
	da[0][1] = 1.;   // differential for a

	// compute differential for f
	double *dfMemory = new double [2 * m];
	double **df      = new double * [m];
	for(i = 0; i < m; i++)
		df[i] = dfMemory + 2 * i;
	forward(tag, int(m), n, d, keep, da, df);

	// print results
	std::cout << "f(a)  = a * (1, 2)^T" << std::endl;
	std::cout << "f'(a) =     (" 
                  << df[0][1] 
                  << ", " 
                  << df[1][1] 
                  << ")^T" 
                  << std::endl;
	bool ok = df[0][1] == 1. && df[1][1] == 2.;

	// free allocated memory
	delete [] daMemory;
	delete [] dfMemory;
	delete [] da;
	delete [] df;


	if( ok ) 
		return 0;    // no error
	return 1;         // error
}

 

// standard input output library
# include <iostream>

// f2cad set up for CppAD::AD<double>
# define F2CAD_USE_CPPAD

// prototype for dscal
# include <f2cad/dscal.hpp>

int main(void)
{	using namespace CppAD;

	// independent variables
	size_t n     = 1;                   // number of independent variables
	CppADvector< AD<double> > a(n);     // vector of independent variables
	a[0]          = 5.;                 // value of independent variables
	Independent(a);                     // declare independent variables

	// create data structure expected by dscal
	AD<double> A[1];
	A[0]         = a[0];
	integer M    = 2;
	AD<double>  F[2];
	F[0]         = 1.;
	F[1]         = 2.;
	integer Incf = 1;

	// set f = a * f 
	f2cad::dscal_(&M, A, F, &Incf );

	// dependent variables
	size_t m = 2;                       // number of dependent variables
	CppADvector< AD<double> > f(2);     // vector of dependent variables
	f[0] = F[0];                        // value of dependent variables
	f[1] = F[1];
	ADFun<double> Fun(a, f);            // declare dependent variables

	// set differential for a 
	CppADvector<double> da(n);
	CppADvector<double> df(m);
	da[0] = 1.;

	// compute differential for f
	df    = Fun.Forward(1, da);

	// print results
	std::cout << "f(a)  = a * (1, 2)^T" << std::endl;
	std::cout << "f'(a) =     (" 
                  << df[0] 
                  << ", " 
                  << df[1] 
                  << ")^T" 
                  << std::endl;

	if( df[0] != 1. || df[1] != 2. )
		return 1;  // error
	return 0;          // no error
}

 

// standard input output library
# include <iostream>

// f2cad set up for F<double>
# define F2CAD_USE_FADBAD

// prototype for dscal
# include <f2cad/dscal.hpp>

int main(void)
{	
	// independent variables
	doublereal  a[1];  // vector of independent variables
	a[0]   = 5.;      // value of independent variables
	a[0].diff(0, 1);  // declare independent variables

	// create data structure expected by dscal
	integer  m = 2;
	doublereal  f[2];
	f[0]     = 1.;
	f[1]     = 2.;
	integer incf = 1;

	// set f = a * f 
	f2cad::dscal_(&m, a, f, &incf );

	// print results
	std::cout << "f(a)  = a * (1, 2)^T" << std::endl;
	std::cout << "f'(a) =     (" 
                  << f[0].d(0) 
                  << ", " 
                  << f[1].d(0) 
                  << ")^T" 
                  << std::endl;

	// return value not significant
	if( f[0].d(0) != 1. || f[1].d(0) != 2. )
		return 1; // error
	return 0;         // no error
}

 
     
     typedef enum { test_none, test_pass, test_fail } test_result;
     
     typedef long int integer;
     
     # ifdef F2CAD_USE_ADOLC
     # include <adolc/adouble.h>
     # include <adolc/interfaces.h>
     # include <adolc/taping.h>
     typedef adouble doublereal;
     # endif
     
     # ifdef F2CAD_USE_CPPAD
     # include <cppad/cppad.hpp>
     typedef CppAD::AD<double> doublereal;
     # endif
     
     # ifdef F2CAD_USE_FADBAD
     # include <FADBAD++/fadiff.h>
     typedef fadbad::F<double> doublereal;
     # endif
     

 

# include <cmath>
# include <f2cad/d9knus.hpp>

test_result d9knus(void)
{	bool ok = true;

	// independent variable
	integer    n = 1;          // number of independent variables
	doublereal x[1];           // vector of independent variables
	x[0]         = .75;        // value of independent variable
	f2cad::Independent(n, x);  // declare x as independent variable vector

	// call the routine
	doublereal xnu = .5;
	doublereal bknu, bknu1;
	integer    iswtch;  
	f2cad::d9knus_(&xnu, x+0, &bknu, &bknu1, &iswtch);

	// dependent variables 
	integer     m = 1;          // number of dependent variables
	doublereal  f[1];           // vector of dependent variables
	f[0] = bknu;                // value f[0] = K_nu (x[0]) * exp( x[0] )
	f2cad::Dependent(m, f);     // declare f as dependent variable vector

	// check function value
	double x0    = f2cad::Value( x[0] );
	double f0    = f2cad::Value( f[0] );
	double pi    = 4. * std::atan(1.);
	double check = 1. / sqrt( 2. * x0 / pi );
	ok          &= f2cad::near_equal(f0, check, 1e-10, 1e-10);

	// Evaluate the derivative of f[0] w.r.t x[0]
	double p = f2cad::Partial<doublereal>(0, 0);
	check    = - check * check * check / pi;
	ok      &= f2cad::near_equal(p, check, 1e-10, 1e-10);

	if( ok ) 
		return test_pass;
	return test_fail;
}

      subroutine d9knus(xnu,x,bknu,bknu1,iswtch)
c***begin prologue  d9knus
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c10b3
c***keywords  bessel function,double precision,fnlib,special function
c***author  fullerton, w., (lanl)
c***purpose  computes bessel functions exp(x)*k-sub-xnu(x) and exp(x)*
c            k-sub-xnu+1(x) for 0. .le. xnu .lt. 1.
c***description
c
c compute bessel functions exp(x) * k-sub-xnu (x)  and
c exp(x) * k-sub-xnu+1 (x) for 0.0 .le. xnu .lt. 1.0 .
c
c series for c0k        on the interval  0.          to  2.50000e-01
c                                        with weighted error   2.16e-32
c                                         log weighted error  31.67
c                               significant figures required  30.86
c                                    decimal places required  32.40
c
c series for znu1       on the interval -7.00000e-01 to  0.
c                                        with weighted error   2.45e-33
c                                         log weighted error  32.61
c                               significant figures required  31.85
c                                    decimal places required  33.26
c***references  (none)
c***routines called  d1mach,dcsevl,dgamma,initds,xerror
c***end prologue  d9knus
      double precision xnu, x, bknu, bknu1, alpha(32), beta(32), a(32),
     1  c0kcs(29), znu1cs(20), alnz, aln2, a0, bknud, bknu0,
     2  b0, c0, euler, expx, p1, p2, p3, qq, result, sqpi2, sqrtx, v,
     3  vlnz, xi, xmu, xnusml, xsml, x2n, x2tov, z, ztov, alnsml,
     4  alnbig
      real alneps
      double precision d1mach, dcsevl, dgamma
      data c0k cs(  1) / +.6018305724 2626108387 5774451803 29 d-1     /
      data c0k cs(  2) / -.1536487143 3017286092 9597559431 24 d+0     /
      data c0k cs(  3) / -.1175117600 8210492040 0682292262 13 d-1     /
      data c0k cs(  4) / -.8524878889 1979509827 0484015509 87 d-3     /
      data c0k cs(  5) / -.6132983876 7496791874 0981769221 11 d-4     /
      data c0k cs(  6) / -.4405228124 5510444562 6798895485 05 d-5     /
      data c0k cs(  7) / -.3163124672 8384488192 9154458921 99 d-6     /
      data c0k cs(  8) / -.2271071938 2899588330 6737717933 96 d-7     /
      data c0k cs(  9) / -.1630564460 8077609552 2746205153 60 d-8     /
      data c0k cs( 10) / -.1170693929 9414776568 7560440431 30 d-9     /
      data c0k cs( 11) / -.8405206378 6464437174 5465934137 92 d-11    /
      data c0k cs( 12) / -.6034667011 8979991487 0960507371 98 d-12    /
      data c0k cs( 13) / -.4332696033 5681371952 0459973669 03 d-13    /
      data c0k cs( 14) / -.3110735803 0203546214 6346977722 37 d-14    /
      data c0k cs( 15) / -.2233407822 6736982254 4861334098 40 d-15    /
      data c0k cs( 16) / -.1603514671 6864226300 6357915286 10 d-16    /
      data c0k cs( 17) / -.1151271736 3666556196 0356977053 05 d-17    /
      data c0k cs( 18) / -.8265759174 6836959105 1694790892 58 d-19    /
      data c0k cs( 19) / -.5934548080 6383948172 3334366959 84 d-20    /
      data c0k cs( 20) / -.4260813819 6467143926 4996130239 76 d-21    /
      data c0k cs( 21) / -.3059126686 4812876299 2636983705 42 d-22    /
      data c0k cs( 22) / -.2196354142 6734575224 9755018155 16 d-23    /
      data c0k cs( 23) / -.1576911326 1495836071 1057506847 60 d-24    /
      data c0k cs( 24) / -.1132171393 5950320948 7577310480 56 d-25    /
      data c0k cs( 25) / -.8128624883 4598404082 7923497144 33 d-27    /
      data c0k cs( 26) / -.5836090089 3453226552 8293493159 49 d-28    /
      data c0k cs( 27) / -.4190124162 3610922519 4523377809 05 d-29    /
      data c0k cs( 28) / -.3008373796 0206435069 5305042128 62 d-30    /
      data c0k cs( 29) / -.2159915206 7808647728 3421680898 32 d-31    /
      data znu1cs(  1) / +.2033067569 9419172967 4444001216 911 d+0    /
      data znu1cs(  2) / +.1400779334 1321977106 2943670790 563 d+0    /
      data znu1cs(  3) / +.7916796961 0016135284 0972241972 320 d-2    /
      data znu1cs(  4) / +.3398011825 3210404535 2930092205 750 d-3    /
      data znu1cs(  5) / +.1174197568 8989336666 4507228352 690 d-4    /
      data znu1cs(  6) / +.3393575706 1226168033 3825865475 121 d-6    /
      data znu1cs(  7) / +.8425941769 7621991019 4629891264 803 d-8    /
      data znu1cs(  8) / +.1833366770 2485008918 4748150900 090 d-9    /
      data znu1cs(  9) / +.3549698447 0441631086 3007064469 557 d-11   /
      data znu1cs( 10) / +.6190324964 6988733220 5244342078 407 d-13   /
      data znu1cs( 11) / +.9819645356 8043942496 0346115456 527 d-15   /
      data znu1cs( 12) / +.1428513143 9649047421 1473563005 985 d-16   /
      data znu1cs( 13) / +.1918949218 8782529896 6162467488 436 d-18   /
      data znu1cs( 14) / +.2394309797 3949891416 2313140597 128 d-20   /
      data znu1cs( 15) / +.2788902468 1534735483 5870465474 995 d-22   /
      data znu1cs( 16) / +.3046066506 3303344258 2845214092 865 d-24   /
      data znu1cs( 17) / +.3131732370 4219181577 1564260932 089 d-26   /
      data znu1cs( 18) / +.3041330989 8785495164 5174908005 034 d-28   /
      data znu1cs( 19) / +.2798403846 3683308434 3185097659 733 d-30   /
      data znu1cs( 20) / +.2446371862 7449759648 5238794922 666 d-32   /
      data euler / 0.5772156649 0153286060 6512090082 40d0 /
      data sqpi2 / +1.253314137 3155002512 0788264240 55 d0      /
      data aln2 / 0.6931471805 5994530941 7232121458 18d0 /
      data ntc0k, ntznu1 / 2*0 /
      data xnusml, xsml, alnsml, alnbig / 4*0.d0 /
      data alneps / 0.0 /
c***first executable statement  d9knus
      if (ntc0k.ne.0) go to 10
      eta = 0.1d0*d1mach(3)
      ntc0k = initds (c0kcs, 29, eta)
      ntznu1 = initds (znu1cs, 20, eta)
c
      xnusml = dsqrt (d1mach(3)/8.d0)
      xsml = 0.1d0*d1mach(3)
      alnsml = dlog (d1mach(1))
      alnbig = dlog (d1mach(2))
      alneps = dlog (0.1d0*d1mach(3))
c
 10   if (xnu.lt.0.d0 .or. xnu.ge.1.d0) call xerror (  'd9knus  xnu must
     1 be ge 0 and lt 1', 33, 1, 2)
      if (x.le.0.) call xerror ( 'd9knus  x must be gt 0', 22, 2, 2)
c
      iswtch = 0
      if (x.gt.2.0d0) go to 50
c
c x is small.  compute k-sub-xnu (x) and the derivative of k-sub-xnu (x)
c then find k-sub-xnu+1 (x).  xnu is reduced to the interval (-.5,+.5)
c then to (0., .5), because k of negative order (-nu) = k of positive
c order (+nu).
c
      v = xnu
      if (xnu.gt.0.5d0) v = 1.0d0 - xnu
c
c carefully find (x/2)**xnu and z**xnu where z = x*x/4.
      alnz = 2.d0 * (dlog(x) - aln2)
c
      if (x.gt.xnu) go to 20
      if (-0.5d0*xnu*alnz-aln2-dlog(xnu).gt.alnbig) call xerror ( 'd9knu
     1s  x so small bessel k-sub-xnu overflows', 45, 3, 2)
c
 20   vlnz = v*alnz
      x2tov = dexp (0.5d0*vlnz)
      ztov = 0.0d0
      if (vlnz.gt.alnsml) ztov = x2tov**2
c
      a0 = 0.5d0*dgamma(1.0d0+v)
      b0 = 0.5d0*dgamma(1.0d0-v)
      c0 = -euler
      if (ztov.gt.0.5d0 .and. v.gt.xnusml) c0 = -0.75d0 +
     1  dcsevl ((8.0d0*v)*v-1.0d0, c0kcs, ntc0k)
c
      if (ztov.le.0.5d0) alpha(1) = (a0-ztov*b0)/v
      if (ztov.gt.0.5d0) alpha(1) = c0 - alnz*(0.75d0 +
     1  dcsevl (vlnz/0.35d0+1.0d0, znu1cs, ntznu1))*b0
      beta(1) = -0.5d0*(a0+ztov*b0)
c
      z = 0.0d0
      if (x.gt.xsml) z = 0.25d0*x*x
      nterms = max1 (2.0, 11.0+(8.*sngl(alnz)-25.19-alneps)
     1  /(4.28-sngl(alnz)))
      do 30 i=2,nterms
        xi = i - 1
        a0 = a0/(xi*(xi-v))
        b0 = b0/(xi*(xi+v))
        alpha(i) = (alpha(i-1)+2.0d0*xi*a0)/(xi*(xi+v))
        beta(i) = (xi-0.5d0*v)*alpha(i) - ztov*b0
 30   continue
c
      bknu = alpha(nterms)
      bknud = beta(nterms)
      do 40 ii=2,nterms
        i = nterms + 1 - ii
        bknu = alpha(i) + bknu*z
        bknud = beta(i) + bknud*z
 40   continue
c
      expx = dexp(x)
      bknu = expx*bknu/x2tov
c
      if (-0.5d0*(xnu+1.d0)*alnz-2.0d0*aln2.gt.alnbig) iswtch = 1
      if (iswtch.eq.1) return
      bknud = expx*bknud*2.0d0/(x2tov*x)
c
      if (xnu.le.0.5d0) bknu1 = v*bknu/x - bknud
      if (xnu.le.0.5d0) return
c
      bknu0 = bknu
      bknu = -v*bknu/x - bknud
      bknu1 = 2.0d0*xnu*bknu/x + bknu0
      return
c
c x is large.  find k-sub-xnu (x) and k-sub-xnu+1 (x) with y. l. luke-s
c rational expansion.
c
 50   sqrtx = dsqrt(x)
      if (x.gt.1.0d0/xsml) go to 90
      an = -0.60 - 1.02/sngl(x)
      bn = -0.27 - 0.53/sngl(x)
      nterms = min0 (32, max1 (3.0, an+bn*alneps))
c
      do 80 inu=1,2
        xmu = 0.d0
        if (inu.eq.1 .and. xnu.gt.xnusml) xmu = (4.0d0*xnu)*xnu
        if (inu.eq.2) xmu = 4.0d0*(dabs(xnu)+1.d0)**2
c
        a(1) = 1.0d0 - xmu
        a(2) = 9.0d0 - xmu
        a(3) = 25.0d0 - xmu
        if (a(2).eq.0.d0) result = sqpi2*(16.d0*x+xmu+7.d0) /
     1    (16.d0*x*sqrtx)
        if (a(2).eq.0.d0) go to 70
c
        alpha(1) = 1.0d0
        alpha(2) = (16.d0*x+a(2))/a(2)
        alpha(3) = ((768.d0*x+48.d0*a(3))*x + a(2)*a(3))/(a(2)*a(3))
c
        beta(1) = 1.0d0
        beta(2) = (16.d0*x+(xmu+7.d0))/a(2)
        beta(3) = ((768.d0*x+48.d0*(xmu+23.d0))*x +
     1    ((xmu+62.d0)*xmu+129.d0))/(a(2)*a(3))
c
        if (nterms.lt.4) go to 65
        do 60 i=4,nterms
          n = i - 1
          x2n = 2*n - 1
c
          a(i) = (x2n+2.d0)**2 - xmu
          qq = 16.d0*x2n/a(i)
          p1 = -x2n*(dble(float(12*n*n-20*n))-a(1))/((x2n-2.d0)*a(i))
     1      - qq*x
          p2 = (dble(float(12*n*n-28*n+8))-a(1))/a(i) - qq*x
          p3 = -x2n*a(i-3)/((x2n-2.d0)*a(i))
c
          alpha(i) = -p1*alpha(i-1) - p2*alpha(i-2) - p3*alpha(i-3)
          beta(i) = -p1*beta(i-1) - p2*beta(i-2) - p3*beta(i-3)
 60     continue
c
 65     result = sqpi2*beta(nterms)/(sqrtx*alpha(nterms))
c
 70     if (inu.eq.1) bknu = result
        if (inu.eq.2) bknu1 = result
 80   continue
      return
c
 90   bknu = sqpi2/sqrtx
      bknu1 = bknu
      return
c
      end

 
# include <f2cad/d9lgmc.hpp>

test_result d9lgmc(void)
{
	// Replace the following line by an example / test of d9lgmc
	// that returns test_pass if the test passes and test_fail otherwise.
	return test_none;
}

      double precision function d9lgmc(x)
c***begin prologue  d9lgmc
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c7e
c***keywords  complete gamma function,correction factor,
c             double precision,gamma function,logarithm,
c             special function
c***author  fullerton, w., (lanl)
c***purpose  computes the  d.p. log gamma correction factor for
c            x .ge. 10. so that dlog(dgamma(x)) = dlog(dsqrt(2*pi)) +
c            (x-5.)*dlog(x) - x + d9lgmc(x)
c***description
c
c compute the log gamma correction factor for x .ge. 10. so that
c dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x)
c
c series for algm       on the interval  0.          to  1.00000e-02
c                                        with weighted error   1.28e-31
c                                         log weighted error  30.89
c                               significant figures required  29.81
c                                    decimal places required  31.48
c***references  (none)
c***routines called  d1mach,dcsevl,initds,xerror
c***end prologue  d9lgmc
      double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach
      data algmcs(  1) / +.1666389480 4518632472 0572965082 2 d+0      /
      data algmcs(  2) / -.1384948176 0675638407 3298605913 5 d-4      /
      data algmcs(  3) / +.9810825646 9247294261 5717154748 7 d-8      /
      data algmcs(  4) / -.1809129475 5724941942 6330626671 9 d-10     /
      data algmcs(  5) / +.6221098041 8926052271 2601554341 6 d-13     /
      data algmcs(  6) / -.3399615005 4177219443 0333059966 6 d-15     /
      data algmcs(  7) / +.2683181998 4826987489 5753884666 6 d-17     /
      data algmcs(  8) / -.2868042435 3346432841 4462239999 9 d-19     /
      data algmcs(  9) / +.3962837061 0464348036 7930666666 6 d-21     /
      data algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23     /
      data algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24     /
      data algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26     /
      data algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27     /
      data algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29     /
      data algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30     /
      data nalgm, xbig, xmax / 0, 2*0.d0 /
c***first executable statement  d9lgmc
      if (nalgm.ne.0) go to 10
      nalgm = initds (algmcs, 15, sngl(d1mach(3)) )
      xbig = 1.0d0/dsqrt(d1mach(3))
      xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1))))
c
 10   if (x.lt.10.d0) call xerror ( 'd9lgmc  x must be ge 10', 23, 1, 2)
      if (x.ge.xmax) go to 20
c
      d9lgmc = 1.d0/(12.d0*x)
      if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs,
     1  nalgm) / x
      return
c
 20   d9lgmc = 0.d0
      call xerror ( 'd9lgmc  x so big d9lgmc underflows', 34, 2, 1)
      return
c
      end

 

# include <f2cad/daxpy.hpp>

test_result daxpy(void)
{	bool ok = true;

	// Input values for daxpy
	doublereal a = 5.;

	integer n    = 2;
	doublereal  x[2];
	x[0]         = 1.;
	x[1]         = 2.;
	integer incx = 1;

	integer m    = 2;
	doublereal  f[2];
	f[0]         = 3.;
	f[1]         = 4.;
	integer incf = 1;

	// declare independent variables
	f2cad::Independent(n, x);

	// set f = a * x + f
	f2cad::daxpy_(&n, &a, x, &incx, f, &incf);

	// declare dependent variables
	f2cad::Dependent(m, f);

	double p;
	p   = f2cad::Partial<doublereal>(0, 0);
	ok &= f2cad::near_equal(p, 5., 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(1, 0);
	ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(0, 1);
	ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(1, 1);
	ok &= f2cad::near_equal(p, 5., 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      subroutine daxpy(n,da,dx,incx,dy,incy)
c***begin prologue  daxpy
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a7
c***keywords  blas,double precision,linear algebra,triad,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p computation y = a*x + y
c***description
c
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       da  double precision scalar multiplier
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c       dy  double precision vector with n elements
c     incy  storage spacing between elements of dy
c
c     --output--
c       dy  double precision result (unchanged if n .le. 0)
c
c     overwrite double precision dy with double precision da*dx + dy.
c     for i = 0 to n-1, replace  dy(ly+i*incy) with da*dx(lx+i*incx) +
c       dy(ly+i*incy), where lx = 1 if incx .ge. 0, else lx = (-incx)*n
c       and ly is defined in a similar way using incy.
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  daxpy
c
      double precision dx(1),dy(1),da
c***first executable statement  daxpy
      if(n.le.0.or.da.eq.0.d0) return
      if(incx.eq.incy) if(incx-1) 5,20,60
    5 continue
c
c        code for nonequal or nonpositive increments.
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop so remaining vector length is a multiple of 4.
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
c
c        code for equal, positive, nonunit increments.
c
   60 continue
      ns = n*incx
          do 70 i=1,ns,incx
          dy(i) = da*dx(i) + dy(i)
   70     continue
      return
      end

 

# include <cmath>
# include <f2cad/dbesks.hpp>

test_result dbesks(void)
{	bool ok = true;

	// independent variable
	integer    n = 1;          // number of independent variables
	doublereal x[1];           // vector of independent variables
	x[0]         = .75;        // value of independent variable
	f2cad::Independent(n, x);  // declare x as independent variable vector

	// call the routine
	doublereal xnu = .5;
	doublereal bk[2];
	integer    nin = 2;  
	f2cad::dbesks_(&xnu, x+0, &nin, bk);

	// dependent variables 
	integer     m = 1;          // number of dependent variables
	doublereal  f[1];           // vector of dependent variables
	f[0] = bk[1];               // value f[0] = K_(nu+1) (x[0])
	f2cad::Dependent(m, f);     // declare f as dependent variable vector

	// check function value
	double x0    = f2cad::Value( x[0] );
	double f0    = f2cad::Value( f[0] );
	double pi    = 4. * std::atan(1.);
	double term  = 1. / sqrt( 2. * x0 / pi );
	double check = term * (1. + 1. / x0 ) * exp(-x0);
	ok          &= f2cad::near_equal(f0, check, 1e-10, 1e-10);

	// Evaluate the derivative of f[0] w.r.t x[0]
	double p = f2cad::Partial<doublereal>(0, 0);
	check    = - check;
	check   += - term * exp(-x0) / (x0 * x0);
	check   += - (1. + 1. / x0) * exp(-x0) * term * term * term / pi;
	ok      &= f2cad::near_equal(p, check, 1e-10, 1e-10);

	if( ok ) 
		return test_pass;
	return test_fail;
}

      subroutine dbesks(xnu,x,nin,bk)
c***begin prologue  dbesks
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c10b3
c***keywords  bessel function,double precision,fractional order,
c             modified bessel function,sequence,special function,
c             third kind
c***author  fullerton, w., (lanl)
c***purpose  computes a d.p. sequence of modified bessel functions of
c            the third kind of fractional order.
c***description
c
c dbesks computes a sequence of modified bessel functions of the third
c kind of order xnu + i at x, where x .gt. 0, xnu lies in (-1,1),
c and i = 0, 1, ... , nin - 1, if nin is positive and i = 0, 1, ... ,
c nin + 1, if nin is negative.  on return, the vector bk(.) contains
c the results at x for order starting at xnu.  xnu, x, and bk are
c double precision.  nin is an integer.
c***references  (none)
c***routines called  d1mach,dbskes,xerror
c***end prologue  dbesks
      double precision xnu, x, bk(1), expxi, xmax, d1mach
      data xmax / 0.d0 /
c***first executable statement  dbesks
      if (xmax.eq.0.d0) xmax = -dlog (d1mach(1))
c
      if (x.gt.xmax) call xerror ( 'dbesks  x so big bessel k underflows
     1', 36, 1, 2)
c
      call dbskes (xnu, x, nin, bk)
c
      expxi = dexp (-x)
      n = iabs (nin)
      do 20 i=1,n
        bk(i) = expxi * bk(i)
 20   continue
c
      return
      end

 

# include <cmath>
# include <f2cad/dbskes.hpp>

test_result dbskes(void)
{	bool ok = true;

	// independent variable
	integer    n = 1;          // number of independent variables
	doublereal x[1];           // vector of independent variables
	x[0]         = .75;        // value of independent variable
	f2cad::Independent(n, x);  // declare x as independent variable vector

	// call the routine
	doublereal xnu = .5;
	doublereal bke[2];
	integer    nin = 2;  
	f2cad::dbskes_(&xnu, x+0, &nin, bke);

	// dependent variables 
	integer     m = 1;          // number of dependent variables
	doublereal  f[1];           // vector of dependent variables
	f[0] = bke[1];              // value f[0] = K_(nu+1) (x[0]) * exp( x[0] )
	f2cad::Dependent(m, f);     // declare f as dependent variable vector

	// check function value
	double x0    = f2cad::Value( x[0] );
	double f0    = f2cad::Value( f[0] );
	double pi    = 4. * std::atan(1.);
	double term  = 1. / sqrt( 2. * x0 / pi );
	double check = term * (1. + 1. / x0 );
	ok          &= f2cad::near_equal(f0, check, 1e-10, 1e-10);

	// Evaluate the derivative of f[0] w.r.t x[0]
	double p = f2cad::Partial<doublereal>(0, 0);
	check    = - (1. + 1. / x0) * term * term * term / pi;
	check   += - term / (x0 * x0);
	ok      &= f2cad::near_equal(p, check, 1e-10, 1e-10);

	if( ok ) 
		return test_pass;
	return test_fail;
}

      subroutine dbskes(xnu,x,nin,bke)
c***begin prologue  dbskes
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c10b3
c***keywords  bessel function,double precision,exponentially scaled,
c             fractional order,modified bessel function,sequence,
c             special function,third kind
c***author  fullerton, w., (lanl)
c***purpose  computes a d.p. sequence of exponentially scaled modified
c            bessel functions of the third kind of fractional order.
c***description
c
c dbskes(xnu,x,nin,bke) computes a double precision sequence
c of exponentially scaled modified bessel functions
c of the third kind of order xnu + i at x, where x .gt. 0,
c xnu lies in (-1,1), and i = 0, 1, ... , nin - 1, if nin is positive
c and i = 0, -1, ... , nin + 1, if nin is negative.  on return, the
c vector bke(.) contains the results at x for order starting at xnu.
c xnu, x, and bke are double precison.  nin is integer.
c***references  (none)
c***routines called  d1mach,d9knus,xerror
c***end prologue  dbskes
      double precision xnu, x, bke(1), bknu1, v, vincr, vend, alnbig,
     1  d1mach
      data alnbig / 0.d0 /
c***first executable statement  dbskes
      if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2))
c
      v = dabs(xnu)
      n = iabs(nin)
c
      if (v.ge.1.d0) call xerror ( 'dbskes  abs(xnu) must be lt 1', 29,
     1  2, 2)
      if (x.le.0.d0) call xerror ( 'dbskes  x is le 0', 17, 3, 2)
      if (n.eq.0) call xerror ( 'dbskes  n the number in the sequence is
     1 0', 41, 4, 2)
c
      call d9knus (v, x, bke(1), bknu1, iswtch)
      if (n.eq.1) return
c
      vincr = sign (1.0, float(nin))
      direct = vincr
      if (xnu.ne.0.d0) direct = vincr*dsign(1.d0, xnu)
      if (iswtch.eq.1 .and. direct.gt.0.) call xerror ( 'dbskes  x so sm
     1all bessel k-sub-xnu+1 overflows', 47, 5, 2)
      bke(2) = bknu1
c
      if (direct.lt.0.) call d9knus (dabs(xnu+vincr), x, bke(2), bknu1,
     1  iswtch)
      if (n.eq.2) return
c
      vend = dabs (xnu+dble(float(nin))) - 1.0d0
      if ((vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)-.694d0).gt.alnbig)
     1  call xerror (  'dbskes  x so small or abs(nu) so big that bessel
     2 k-sub-nu overflows', 67, 5, 2)
c
      v = xnu
      do 10 i=3,n
        v = v + vincr
        bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2)
 10   continue
c
      return
      end

 

# include <f2cad/dcsevl.hpp>

test_result dcsevl(void)
{	bool ok = true;
	// independent variable
	integer  N = 1;
	double   x = .5;
	doublereal X[1];
	X[0]       = x;
	// dependent variable
	integer  M = 1;
	doublereal F[1];
	// polynomial coefficients
	integer n  = 3;
	double     a[3];
	doublereal A[3];
	for(integer i = 0; i < n; i++)
	{	a[i] = 1. / double(i+1);
		A[i] = a[i];
	}

	// declare independent variables
	f2cad::Independent(N, X);

	F[0] = f2cad::dcsevl_(X, A, &n);

	// declare dependent variables
	f2cad::Dependent(M, F);

	// evaluate the Chebyshev polynomials
	double t_0 = 1.;
	double t_1 = x;
	double t_2 = 2. * x * t_1 - t_0;

	// check function value
	double value = f2cad::Value( F[0] );
	double check = a[0] * t_0 / 2. + a[1] * t_1 + a[2] * t_2;
	ok          &= f2cad::near_equal(check, value, 1e-10, 1e-10);

	// evaluate the Derivatives
	double dt_0 = 0.;
	double dt_1 = 1.;
	double dt_2 = 2. * t_1 + 2. * x * dt_1 - dt_0; 
	value       = f2cad::Partial<doublereal>(0, 0);
	check       = a[1] * dt_1 + a[2] * dt_2;
	ok          &= f2cad::near_equal(check, value, 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      double precision function dcsevl(x,a,n)
c***begin prologue  dcsevl
c***date written   770401   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c3a2
c***keywords  chebyshev,fnlib,special function
c***author  fullerton, w., (lanl)
c***purpose  evaluate the double precision n-term chebyshev series a
c            at x.
c***description
c
c evaluate the n-term chebyshev series a at x.  adapted from
c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).
c w. fullerton, c-3, los alamos scientific laboratory.
c
c       input arguments --
c x    double precision value at which the series is to be evaluated.
c a    double precision array of n terms of a chebyshev series.  in
c      evaluating a, only half of the first coefficient is summed.
c n    number of terms in array a.
c***references  (none)
c***routines called  xerror
c***end prologue  dcsevl
c
       double precision a(n),x,twox,b0,b1,b2
c***first executable statement  dcsevl
       if(n.lt.1)call xerror( 'dcsevl  number of terms le 0', 28, 2,2)
       if(n.gt.1000) call xerror ( 'dcsevl  number of terms gt 1000',
     1   31, 3, 2)
       if ((x.lt.-1.d0) .or. (x.gt.1.d0)) call xerror ( 'dcsevl  x outsi
     1de (-1,+1)', 25, 1, 1)
c
       twox = 2.0d0*x
       b1 = 0.d0
       b0=0.d0
       do 10 i=1,n
         b2=b1
         b1=b0
         ni = n - i + 1
         b0 = twox*b1 - b2 + a(ni)
 10    continue
c
       dcsevl = 0.5d0 * (b0-b2)
c
       return
      end

 

# include <f2cad/ddot.hpp>

test_result ddot(void)
{	bool ok = true;

	integer n   = 4;
	doublereal x[4];
	x[0]        = 1.;
	x[1]        = 2.;
	x[2]        = 3.;
	x[3]        = 4.;

	// declare independent variables
	f2cad::Independent(n, x);

	// set f = x[0] * x[2] + x[1] * x[3] 
	integer m    = 1;
	doublereal f[1];
	integer incx = 1;
	integer n2   = n / 2;
	f[0]         = f2cad::ddot_ (&n2, x, &incx, x+n2, &incx);

	// declare dependent variables
	f2cad::Dependent(m, f);

	double p;
	p   = f2cad::Partial<doublereal>(0, 0);
	ok &= f2cad::near_equal(p, f2cad::Value(x[2]), 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(0, 1);
	ok &= f2cad::near_equal(p, f2cad::Value(x[3]), 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(0, 2);
	ok &= f2cad::near_equal(p, f2cad::Value(x[0]), 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(0, 3);
	ok &= f2cad::near_equal(p, f2cad::Value(x[1]), 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      double precision function ddot(n,dx,incx,dy,incy)
c***begin prologue  ddot
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a4
c***keywords  blas,double precision,inner product,linear algebra,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p. inner product of d.p. vectors
c***description
c
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c       dy  double precision vector with n elements
c     incy  storage spacing between elements of dy
c
c     --output--
c     ddot  double precision dot product (zero if n .le. 0)
c
c     returns the dot product of double precision dx and dy.
c     ddot = sum for i = 0 to n-1 of  dx(lx+i*incx) * dy(ly+i*incy)
c     where lx = 1 if incx .ge. 0, else lx = (-incx)*n, and ly is
c     defined in a similar way using incy.
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  ddot
c
      double precision dx(1),dy(1)
c***first executable statement  ddot
      ddot = 0.d0
      if(n.le.0)return
      if(incx.eq.incy) if(incx-1) 5,20,60
    5 continue
c
c         code for unequal or nonpositive increments.
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
         ddot = ddot + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1.
c
c
c        clean-up loop so remaining vector length is a multiple of 5.
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         ddot = ddot + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
         ddot = ddot + dx(i)*dy(i) + dx(i+1)*dy(i+1) +
     1   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
      return
c
c         code for positive equal increments .ne.1.
c
   60 continue
      ns = n*incx
          do 70 i=1,ns,incx
          ddot = ddot + dx(i)*dy(i)
   70     continue
      return
      end

 

# include <f2cad/dgamlm.hpp>

test_result dgamlm(void)
{	bool ok = true;
	doublereal ad_xmin, ad_xmax;

	f2cad::dgamlm_(&ad_xmin, &ad_xmax);

	double xmin = f2cad::Value(ad_xmin);
	double xmax = f2cad::Value(ad_xmax);

	ok  &= xmin < 0.;
	ok  &= xmax > 0.;
	ok  &= (xmax - xmin) > 50.;
	ok  &= f2cad::near_equal( log(exp(xmin)), xmin, 1e-10, 1e-10); 
	ok  &= f2cad::near_equal( log(exp(xmax)), xmax, 1e-10, 1e-10); 

	if( ok )
		return test_pass;
	return test_fail;
}

      subroutine dgamlm(xmin,xmax)
c***begin prologue  dgamlm
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c7a,r2
c***keywords  complete gamma function,double precision,gamma function,
c             limits,special function
c***author  fullerton, w., (lanl)
c***purpose  computes the d.p. minimum and maximum bounds for x in
c            gamma(x).
c***description
c
c calculate the minimum and maximum legal bounds for x in gamma(x).
c xmin and xmax are not the only bounds, but they are the only non-
c trivial ones to calculate.
c
c             output arguments --
c xmin   double precision minimum legal value of x in gamma(x).  any
c        smaller value of x might result in underflow.
c xmax   double precision maximum legal value of x in gamma(x).  any
c        larger value of x might cause overflow.
c***references  (none)
c***routines called  d1mach,xerror
c***end prologue  dgamlm
      double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach
c***first executable statement  dgamlm
      alnsml = dlog(d1mach(1))
      xmin = -alnsml
      do 10 i=1,10
        xold = xmin
        xln = dlog(xmin)
        xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml)
     1    / (xmin*xln+0.5d0)
        if (dabs(xmin-xold).lt.0.005d0) go to 20
 10   continue
      call xerror ( 'dgamlm  unable to find xmin', 27, 1, 2)
c
 20   xmin = -xmin + 0.01d0
c
      alnbig = dlog (d1mach(2))
      xmax = alnbig
      do 30 i=1,10
        xold = xmax
        xln = dlog(xmax)
        xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig)
     1    / (xmax*xln-0.5d0)
        if (dabs(xmax-xold).lt.0.005d0) go to 40
 30   continue
      call xerror ( 'dgamlm  unable to find xmax', 27, 2, 2)
c
 40   xmax = xmax - 0.01d0
      xmin = dmax1 (xmin, -xmax+1.d0)
c
      return
      end

 
# include <f2cad/dgamma.hpp>

test_result dgamma(void)
{	bool ok = true;

	// independent variables 
	integer     n  = 1;         // number of independent variables
	double     ell = 3.;
	doublereal  x[1];           // vector of independent variables
	x[0]          = ell+1.;     // value  of independent variables
	f2cad::Independent(n, x);   // declare x as independent variable vector

	// dependent variables 
	integer     m = 1;          // number of dependent variables
	doublereal  f[1];           // vector of dependent variables
	f[0] = f2cad::dgamma_(x+0); // compute value f[0] = Gamma(x[0])
	f2cad::Dependent(m, f);     // declare f as dependent variable vector

	// check the value of f
	ok &= f2cad::near_equal(f2cad::Value(f[0]), 6., 1e-10, 1e-10);

	// Compute derivative of f w.r.t x using f2cad common interface.
	// This general interface works with adolc, cppad, and fadbad.
	// See GetStarted for package specific more efficient code.
	double gamma = 0.577215664901532;
	double p     = f2cad::Partial<doublereal>(0, 0);
	double check = 6. * ( - gamma + 1. + 1./2. + 1./3. ); 
	ok &= f2cad::near_equal(p, check, 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      double precision function dgamma(x)
c***begin prologue  dgamma
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c7a
c***keywords  complete gamma function,double precision,gamma function,
c             special function
c***author  fullerton, w., (lanl)
c***purpose  computes the d.p. complete gamma function.
c***description
c
c dgamma(x) calculates the double precision complete gamma function
c for double precision argument x.
c
c series for gam        on the interval  0.          to  1.00000e+00
c                                        with weighted error   5.79e-32
c                                         log weighted error  31.24
c                               significant figures required  30.00
c                                    decimal places required  32.05
c***references  (none)
c***routines called  d1mach,d9lgmc,dcsevl,dgamlm,dint,initds,xerror
c***end prologue  dgamma
      double precision x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax,
     1  xmin, y,  dint, d9lgmc, dcsevl, d1mach
c
      data gam cs(  1) / +.8571195590 9893314219 2006239994 2 d-2      /
      data gam cs(  2) / +.4415381324 8410067571 9131577165 2 d-2      /
      data gam cs(  3) / +.5685043681 5993633786 3266458878 9 d-1      /
      data gam cs(  4) / -.4219835396 4185605010 1250018662 4 d-2      /
      data gam cs(  5) / +.1326808181 2124602205 8400679635 2 d-2      /
      data gam cs(  6) / -.1893024529 7988804325 2394702388 6 d-3      /
      data gam cs(  7) / +.3606925327 4412452565 7808221722 5 d-4      /
      data gam cs(  8) / -.6056761904 4608642184 8554829036 5 d-5      /
      data gam cs(  9) / +.1055829546 3022833447 3182350909 3 d-5      /
      data gam cs( 10) / -.1811967365 5423840482 9185589116 6 d-6      /
      data gam cs( 11) / +.3117724964 7153222777 9025459316 9 d-7      /
      data gam cs( 12) / -.5354219639 0196871408 7408102434 7 d-8      /
      data gam cs( 13) / +.9193275519 8595889468 8778682594 0 d-9      /
      data gam cs( 14) / -.1577941280 2883397617 6742327395 3 d-9      /
      data gam cs( 15) / +.2707980622 9349545432 6654043308 9 d-10     /
      data gam cs( 16) / -.4646818653 8257301440 8166105893 3 d-11     /
      data gam cs( 17) / +.7973350192 0074196564 6076717535 9 d-12     /
      data gam cs( 18) / -.1368078209 8309160257 9949917230 9 d-12     /
      data gam cs( 19) / +.2347319486 5638006572 3347177168 8 d-13     /
      data gam cs( 20) / -.4027432614 9490669327 6657053469 9 d-14     /
      data gam cs( 21) / +.6910051747 3721009121 3833697525 7 d-15     /
      data gam cs( 22) / -.1185584500 2219929070 5238712619 2 d-15     /
      data gam cs( 23) / +.2034148542 4963739552 0102605193 2 d-16     /
      data gam cs( 24) / -.3490054341 7174058492 7401294910 8 d-17     /
      data gam cs( 25) / +.5987993856 4853055671 3505106602 6 d-18     /
      data gam cs( 26) / -.1027378057 8722280744 9006977843 1 d-18     /
      data gam cs( 27) / +.1762702816 0605298249 4275966074 8 d-19     /
      data gam cs( 28) / -.3024320653 7353062609 5877211204 2 d-20     /
      data gam cs( 29) / +.5188914660 2183978397 1783355050 6 d-21     /
      data gam cs( 30) / -.8902770842 4565766924 4925160106 6 d-22     /
      data gam cs( 31) / +.1527474068 4933426022 7459689130 6 d-22     /
      data gam cs( 32) / -.2620731256 1873629002 5732833279 9 d-23     /
      data gam cs( 33) / +.4496464047 8305386703 3104657066 6 d-24     /
      data gam cs( 34) / -.7714712731 3368779117 0390152533 3 d-25     /
      data gam cs( 35) / +.1323635453 1260440364 8657271466 6 d-25     /
      data gam cs( 36) / -.2270999412 9429288167 0231381333 3 d-26     /
      data gam cs( 37) / +.3896418998 0039914493 2081663999 9 d-27     /
      data gam cs( 38) / -.6685198115 1259533277 9212799999 9 d-28     /
      data gam cs( 39) / +.1146998663 1400243843 4761386666 6 d-28     /
      data gam cs( 40) / -.1967938586 3451346772 9510399999 9 d-29     /
      data gam cs( 41) / +.3376448816 5853380903 3489066666 6 d-30     /
      data gam cs( 42) / -.5793070335 7821357846 2549333333 3 d-31     /
      data pi / 3.1415926535 8979323846 2643383279 50 d0 /
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
      data ngam, xmin, xmax, dxrel / 0, 3*0.d0 /
c***first executable statement  dgamma
      if (ngam.ne.0) go to 10
      ngam = initds (gamcs, 42, 0.1*sngl(d1mach(3)) )
c
      call dgamlm (xmin, xmax)
      dxrel = dsqrt (d1mach(4))
c
 10   y = dabs(x)
      if (y.gt.10.d0) go to 50
c
c compute gamma(x) for -xbnd .le. x .le. xbnd.  reduce interval and find
c gamma(1+y) for 0.0 .le. y .lt. 1.0 first of all.
c
      n = x
      if (x.lt.0.d0) n = n - 1
      y = x - dble(float(n))
      n = n - 1
      dgamma = 0.9375d0 + dcsevl (2.d0*y-1.d0, gamcs, ngam)
      if (n.eq.0) return
c
      if (n.gt.0) go to 30
c
c compute gamma(x) for x .lt. 1.0
c
      n = -n
      if (x.eq.0.d0) call xerror ( 'dgamma  x is 0', 14, 4, 2)
      if (x.lt.0.0 .and. x+dble(float(n-2)).eq.0.d0) call xerror ( 'dgam
     1ma  x is a negative integer', 31, 4, 2)
      if (x.lt.(-0.5d0) .and. dabs((x-dint(x-0.5d0))/x).lt.dxrel) call
     1  xerror ( 'dgamma  answer lt half precision because x too near ne
     2gative integer', 68, 1, 1)
c
      do 20 i=1,n
        dgamma = dgamma/(x+dble(float(i-1)) )
 20   continue
      return
c
c gamma(x) for x .ge. 2.0 and x .le. 10.0
c
 30   do 40 i=1,n
        dgamma = (y+dble(float(i))) * dgamma
 40   continue
      return
c
c gamma(x) for dabs(x) .gt. 10.0.  recall y = dabs(x).
c
 50   if (x.gt.xmax) call xerror ( 'dgamma  x so big gamma overflows',
     1  32, 3, 2)
c
      dgamma = 0.d0
      if (x.lt.xmin) call xerror ( 'dgamma  x so small gamma underflows'
     1  , 35, 2, 1)
      if (x.lt.xmin) return
c
      dgamma = dexp ((y-0.5d0)*dlog(y) - y + sq2pil + d9lgmc(y) )
      if (x.gt.0.d0) return
c
      if (dabs((x-dint(x-0.5d0))/x).lt.dxrel) call xerror ( 'dgamma  ans
     1wer lt half precision, x too near negative integer'  , 61, 1, 1)
c
      sinpiy = dsin (pi*y)
      if (sinpiy.eq.0.d0) call xerror ( 'dgamma  x is a negative integer
     1', 31, 4, 2)
c
      dgamma = -pi/(y*sinpiy*dgamma)
c
      return
      end

 

# include <f2cad/dgefa.hpp>
# include <f2cad/dgesl.hpp>

test_result dgefa(void)
{	bool ok = true;

	// A is a 2 by 2 matrix in column major order
	doublereal A[4];
	A[0] = doublereal(1); A[2] = doublereal(2);
	A[1] = doublereal(0); A[3] = doublereal(3); 

	// b is a column vector of length 2
	doublereal b[2];
	b[0] = doublereal(1);
	b[1] = doublereal(2);

	// declare independent variables
	integer n = 2;
	f2cad::Independent(n, b);

	// f is a column vector of length 2
	doublereal f[2];
	integer i;
	for(i = 0; i < n; i++)
		f[i] = b[i];

	// solve A * f = b
	integer ipvt[2];
	integer lda = n;
	integer job = 0;
	integer info;
	f2cad::dgefa_ (A, &lda, &n, ipvt, &info);
	f2cad::dgesl_ (A, &lda, &n, ipvt, f, &job);

	// check info flag
	ok &= (info == 0);

	// construct the AD function object F : b -> f
	// f[1] = b[1] / 3
	// f[0] = b[0] - 2 * b[1] / 3 
	integer m = 2;
	f2cad::Dependent(m, f);

	double p;
	p    = f2cad::Partial<doublereal>(0, 0);
	ok   &= f2cad::near_equal(p, 1., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(0, 1);
	ok   &= f2cad::near_equal(p, -2./3., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(1, 0);
	ok   &= f2cad::near_equal(p, 0., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(1, 1);
	ok   &= f2cad::near_equal(p,  1./3., 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      subroutine dgefa(a,lda,n,ipvt,info)
c***begin prologue  dgefa
c***date written   780814   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d2a1
c***keywords  double precision,factor,linear algebra,linpack,matrix
c***author  moler, c. b., (u. of new mexico)
c***purpose  factors a double precision matrix by gaussian elimination.
c***description
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c***references  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c                 *linpack users  guide*, siam, 1979.
c***routines called  daxpy,dscal,idamax
c***end prologue  dgefa
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
c
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c     gaussian elimination with partial pivoting
c
c***first executable statement  dgefa
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end

 

# include <f2cad/dgefa.hpp>
# include <f2cad/dgesl.hpp>

test_result dgesl(void)
{	bool ok = true;

	// A is a 2 by 2 matrix in column major order
	doublereal A[4];
	A[0] = doublereal(1); A[2] = doublereal(2);
	A[1] = doublereal(3); A[3] = doublereal(4); 

	// b is a column vector of length 2
	doublereal b[2];
	b[0] = doublereal(1);
	b[1] = doublereal(2);

	// declare independent variables
	integer n = 2;
	f2cad::Independent(n, b);

	// f is a column vector of length 2
	doublereal f[2];
	integer i;
	for(i = 0; i < n; i++)
		f[i] = b[i];

	// solve A * f = b
	integer ipvt[2];
	integer lda = n;
	integer job = 0;
	integer info;
	f2cad::dgefa_ (A, &lda, &n, ipvt, &info);
	f2cad::dgesl_ (A, &lda, &n, ipvt, f, &job);

	// check info flag
	ok &= (info == 0);

	// construct the AD function object F : b -> f
	// f[0] = b[1] - 2 b[0]
	// f[1] = 3 b[0] / 2 - b[1]/ 2
	integer m = 2;
	f2cad::Dependent(m, f);

	double p;
	p    = f2cad::Partial<doublereal>(0, 0);
	ok   &= f2cad::near_equal(p, -2., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(0, 1);
	ok   &= f2cad::near_equal(p, 1., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(1, 0);
	ok   &= f2cad::near_equal(p, 3./2., 1e-10, 1e-10);

	p    = f2cad::Partial<doublereal>(1, 1);
	ok   &= f2cad::near_equal(p,  -1./2., 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      subroutine dgesl(a,lda,n,ipvt,b,job)
c***begin prologue  dgesl
c***date written   780814   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d2a1
c***keywords  double precision,linear algebra,linpack,matrix,solve
c***author  moler, c. b., (u. of new mexico)
c***purpose  solves the double precision system  a*x=b or  trans(a)*x=b
c            using the factors computed by dgeco or dgefa.
c***description
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c***references  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c                 *linpack users  guide*, siam, 1979.
c***routines called  daxpy,ddot
c***end prologue  dgesl
      integer lda,n,ipvt(1),job
      double precision a(lda,1),b(1)
c
      double precision ddot,t
      integer k,kb,l,nm1
c***first executable statement  dgesl
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end

 
# include <f2cad/dint.hpp>

test_result dint(void)
{	bool ok = true;

	// independent variables 
	integer     n = 2;        // number of independent variables
	doublereal  x[2];         // vector of independent variables
	x[0]          =  4.5;     // value  of independent variables
	x[1]          = -3.5;
	f2cad::Independent(n, x); // declare x as independent variable vector

	// dependent variables 
	integer     m = n;           // number of dependent variables
	doublereal  f[2];            // vector of dependent variables
	integer i;
	for(i = 0; i < m; i++)
		f[i] = f2cad::dint_(x+i);
	f2cad::Dependent(m, f);      // declare f as dependent variable vector

	// check function value
	ok &= (f2cad::Value(f[0]) ==  4.);
	ok &= (f2cad::Value(f[1]) == -3.);

	// Compute derivative of f w.r.t x using f2cad common interface.
	// This general interface works with adolc, cppad, and fadbad.
	// See GetStarted for package specific more efficient code.
	for(i = 0; i < m; i++)
	{	double p = f2cad::Partial<doublereal>(i, i);
		ok      &= f2cad::near_equal(p, 0., 1e-10, 1e-10);
	}

	if( ok )
		return test_pass;
	return test_fail;
}

       double precision function dint (x)
c***begin prologue   dint
c***revision  october 1, 1980
c***category no.  m2
c***keyword(s) truncation,greatest integer,double precision
c***author  fullerton w. (lasl)
c***date written  august 1979
c***purpose  to find the largest integer whose magnitude does
c            not exceed x and convert to double precision.
c***description
c august 1979 edition. w. fullerton, los alamos scientific laboratory.
c installed on the vax by dolores montano, c-3, 5/80.
c
c dint is the double precision equivalent of aint. this portable
c version is quite efficient when the argument is reasonably small (a
c common case), and so no faster machine-dependent version is needed.
c
c
c***reference(s)
c***routines called  d1mach,i1mach,r1mach,xerror
c***end prologue
       double precision x,xscl,scale,xbig,xmax,part,d1mach
       data npart,scale,xbig,xmax /0, 3*0.0d0 /
c***first executable statement    dint
       if(npart.ne.0) go to 10
       ibase = i1mach(7)
       ndigd = i1mach(14)
       ndigi = min0(i1mach(8),i1mach(11)) - 1
       nmax  = 1.0d0/d1mach(4)
       xbig  = amin1 (float(i1mach(9)),1.0/r1mach(4))
       if (ibase.ne.i1mach(10)) call xerror (
     1   'dint   algorithm error. integer base ne real base', 49, 2, 2)
c
       npart = (ndigd + ndigi -1) / ndigi
       scale = ibase**ndigi
c
 10    if(x.lt.(-xbig) .or. x.gt.xbig) go to 20
c
       dint=int(sngl(x))
       return
c
 20    xscl = dabs(x)
       if(xscl.gt.xmax) go to 50
c
       do 30 i=1,npart
       xscl=xscl/scale
 30    continue
c
       do 40 i=1,npart
       xscl = xscl*scale
       ipart=xscl
       part=ipart
       xscl=xscl - part
       dint = dint * scale + part
 40    continue
c
       if(x.lt.0.0d0) dint = -dint
       return
c
 50     call xerror ('dint    dabs(x) may be too big to be represented
     1as an exact integer', 67, 1, 1)
       dint=x
       return
c
       end

 

# include <f2cad/dscal.hpp>

test_result dscal(void)
{	bool ok = true;

	// Input values for dscal
	integer n     = 1;
	doublereal  a[1];
	a[0]          = 5.;

	// declare independent variables
	f2cad::Independent(n, a);

	integer m    = 2;
	doublereal  f[2];
	f[0]         = 1.;
	f[1]         = 2.;
	integer incf = 1;

	// set f = a * f 
	f2cad::dscal_(&n, a, f, &incf );

	// declare dependent variables
	f2cad::Dependent(m, f);

	double p;
	p   = f2cad::Partial<doublereal>(0, 0);
	ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10);

	p   = f2cad::Partial<doublereal>(1, 0);
	ok &= f2cad::near_equal(p, 2., 1e-10, 1e-10);

	if( ok )
		return test_pass;
	return test_fail;
}

      subroutine dscal(n,da,dx,incx)
c***begin prologue  dscal
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a6
c***keywords  blas,linear algebra,scale,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p. vector scale x = a*x
c***description
c
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       da  double precision scale factor
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c
c     --output--
c       dx  double precision result (unchanged if n.le.0)
c
c     replace double precision dx by double precision da*dx.
c     for i = 0 to n-1, replace dx(1+i*incx) with  da * dx(1+i*incx)
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  dscal
c
      double precision da,dx(1)
c***first executable statement  dscal
      if(n.le.0)return
      if(incx.eq.1)goto 20
c
c        code for increments not equal to 1.
c
      ns = n*incx
          do 10 i = 1,ns,incx
          dx(i) = da*dx(i)
   10     continue
      return
c
c        code for increments equal to 1.
c
c
c        clean-up loop so remaining vector length is a multiple of 5.
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end

 

# include <f2cad/idamax.hpp>

test_result idamax(void)
{	bool ok = true;

	// dx is a vector of length 4
	double     data[] = {1, 2, 5, 4};
	doublereal dx[4];
	integer i;
	for(i = 0; i < 4; i++)
		dx[i] = data[i];

	// other arguments to idamax
	integer n    = 4;
	integer incx = 1;

	// check return value
	ok &= ( 3 == f2cad::idamax_(&n, dx, &incx) );

	if( ok )
		return test_pass;
	return test_fail;
}

      integer function idamax(n,dx,incx)
c***begin prologue  idamax
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a2
c***keywords  blas,double precision,linear algebra,maximum component,
c             vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  find largest component of d.p. vector
c***description
c
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c
c     --output--
c   idamax  smallest index (zero if n .le. 0)
c
c     find smallest index of maximum magnitude of double precision dx.
c     idamax =  first i, i = 1 to n, to minimize  abs(dx(1-incx+i*incx)
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  idamax
c
      double precision dx(1),dmax,xmag
c***first executable statement  idamax
      idamax = 0
      if(n.le.0) return
      idamax = 1
      if(n.le.1)return
      if(incx.eq.1)goto 20
c
c        code for increments not equal to 1.
c
      dmax = dabs(dx(1))
      ns = n*incx
      ii = 1
          do 10 i = 1,ns,incx
          xmag = dabs(dx(i))
          if(xmag.le.dmax) go to 5
          idamax = ii
          dmax = xmag
    5     ii = ii + 1
   10     continue
      return
c
c        code for increments equal to 1.
c
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
          xmag = dabs(dx(i))
          if(xmag.le.dmax) go to 30
          idamax = i
          dmax = xmag
   30 continue
      return
      end

 

# include <f2cad/initds.hpp>

test_result initds(void)
{	bool ok = true;

	//                   a_0,  a_1,  a_2,  a_3,  a_4,  a_5
	doublereal a[]   = { 1e0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5};
	integer    n     = sizeof(a) / sizeof(a[0]);
	doublereal eta   = 1e-3;
	integer    m     = f2cad::initds_(a, &n, &eta);
	
	ok &= (m == 4);

	if( ok )
		return test_pass;
	return test_fail;
}

      function initds(dos,nos,eta)
c***begin prologue  initds
c***date written   770601   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  c3a2
c***keywords  chebyshev,double precision,initialize,
c             orthogonal polynomial,series,special function
c***author  fullerton, w., (lanl)
c***purpose  initializes the d.p. properly normalized orthogonal
c            polynomial series to determine the number of terms needed
c            for specific accuracy.
c***description
c
c initialize the double precision orthogonal series dos so that initds
c is the number of terms needed to insure the error is no larger than
c eta.  ordinarily eta will be chosen to be one-tenth machine precision
c
c             input arguments --
c dos    dble prec array of nos coefficients in an orthogonal series.
c nos    number of coefficients in dos.
c eta    requested accuracy of series.
c***references  (none)
c***routines called  xerror
c***end prologue  initds
c
      double precision dos(nos)
c***first executable statement  initds
      if (nos.lt.1) call xerror ( 'initds  number of coefficients lt 1',
     1 35, 2, 2)
c
      err = 0.
      do 10 ii=1,nos
        i = nos + 1 - ii
        err = err + abs(sngl(dos(i)))
        if (err.gt.eta) go to 20
 10   continue
c
 20   if (i.eq.nos) call xerror ( 'initds  eta may be too small', 28,
     1  1, 2)
      initds = i
c
      return
      end

 
	Ok:    xerror

 

# include <cstring>
# include <f2cad/xerror.hpp>

void xerror(void)
{
	const char* messg     = "Pass: xerror\n";
	integer     nmessg    = std::strlen(messg);
	integer     nerr      = 1;
	integer     level     = 0;
 	ftnlen      messg_len = (ftnlen) nmessg;   

	f2cad::xerror_(messg, &nmessg, &nerr, &level, messg_len);

}

      subroutine xerror(messg,nmessg,nerr,level)
c***begin prologue  xerror
c***date written   790801   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  r3c
c***keywords  error,xerror package
c***author  jones, r. e., (snla)
c***purpose  processes an error (diagnostic) message.
c***description
c     abstract
c        xerror processes a diagnostic message, in a manner
c        determined by the value of level and the current value
c        of the library error control flag, kontrl.
c        (see subroutine xsetf for details.)
c
c     description of parameters
c      --input--
c        messg - the hollerith message to be processed, containing
c                no more than 72 characters.
c        nmessg- the actual number of characters in messg.
c        nerr  - the error number associated with this message.
c                nerr must not be zero.
c        level - error category.
c                =2 means this is an unconditionally fatal error.
c                =1 means this is a recoverable error.  (i.e., it is
c                   non-fatal if xsetf has been appropriately called.)
c                =0 means this is a warning message only.
c                =-1 means this is a warning message which is to be
c                   printed at most once, regardless of how many
c                   times this call is executed.
c
c     examples
c        call xerror('smooth -- num was zero.',23,1,2)
c        call xerror('integ  -- less than full accuracy achieved.',
c                    43,2,1)
c        call xerror('rooter -- actual zero of f found before interval f
c    1ully collapsed.',65,3,0)
c        call xerror('exp    -- underflows being set to zero.',39,1,-1)
c
c     latest revision ---  19 mar 1980
c     written by ron jones, with slatec common math library subcommittee
c***references  jones r.e., kahaner d.k., "xerror, the slatec error-
c                 handling package", sand82-0800, sandia laboratories,
c                 1982.
c***routines called  xerrwv
c***end prologue  xerror
      character*(*) messg
c***first executable statement  xerror
c     call xerrwv(messg,nmessg,nerr,level,0,0,0,0,0.,0.)
c statement above was commented out and the modified implementation
c below is for standard linking to a c routine
c
c -------------------------------------------------------------------------
c This is a limited version of the specifications above for use with f2cad
      integer message(1000)
      n = min(nmessg, 1000)
      do i = 1, n
          message(i) = ichar( messg(i:i) )
100   enddo
      call xerror2cpp(message, nmessg, nerr, level) 
      return
      end

 
	All the tests passed.

 
	At least one test failed.

 

// system include files used for I/O
# include <iostream>

// define the types integer and doublereal
# include <f2cad/f2cad.hpp>

// ../add2lib.sh uses this comment for extern command 
extern test_result dbesks(void);
extern test_result dbskes(void);
extern test_result d9knus(void);
extern test_result d9lgmc(void);
extern test_result dgamma(void);
extern test_result dint(void);
extern test_result dcsevl(void);
extern test_result dgamlm(void);
extern test_result initds(void);
extern test_result idamax(void);
extern test_result daxpy(void);
extern test_result ddot(void);
extern test_result dgefa(void);
extern test_result dgesl(void);
extern test_result dscal(void);

// xerror is a special case
extern void xerror(void);

// function that runs one test
bool run(test_result test_fun(void), const char *name)
{
	bool ok = true;
	switch( test_fun() )
	{
		case test_none:
		std::cout << "None: " << name << std::endl;
		break;

		case test_pass:
		std::cout << "Pass: " << name << std::endl;
		break;

		case test_fail:
		std::cout << "Fail: " << name << std::endl;
		ok = false;
		break;
	}

	return ok;
}

// main program that runs all the examples
int main(void)
{	bool ok = true;
	using std::cout;
	using std::endl;

	// ../add2lib.sh uses this comment for run command 
	ok &= run( dbesks , "dbesks" );
	ok &= run( dbskes , "dbskes" );
	ok &= run( d9knus , "d9knus" );
	ok &= run( d9lgmc , "d9lgmc" );
	ok &= run( dgamma , "dgamma" );
	ok &= run( dint , "dint" );
	ok &= run( dcsevl , "dcsevl" );
	ok &= run( dgamlm , "dgamlm" );
	ok &= run( initds , "initds" );
	ok &= run( idamax , "idamax" );
	ok &= run( daxpy  , "daxpy"  );
	ok &= run( ddot   , "ddot"   );
	ok &= run( dgefa  , "dgefa"  );
	ok &= run( dgesl  , "dgesl"  );
	ok &= run( daxpy  , "dscal"  );

	if( ok )
		cout << endl << "Pass: None of the tests above failed." << endl;
	else	cout << endl << "Fail: At least one test above failed." << endl;
	cout << "      The next two lines should be identical:" << endl;

	// xerror is a special case because it is not checked automatically
	// and must be checked by eye.
	cout << "Pass: xerror" << endl;
	xerror();
	cout << endl;

	if( ok )
		return 0;
	return 1;
}

 // assume that F2CAD_USE_Adolc is defined
# include <f2cad/f2cad.hpp>

namespace {
	static bool Recording  = false; // currently recording a function
	static int Domain      = 0;     // domain dimension
	static int Range       = 0;     // range dimension
	static double *Xmemory = 0;     // domain value and differential
	static double **X      = 0;
	static double *Ymemory = 0;     // range value and differential
	static double **Y      = 0;
}

namespace f2cad {

void Independent(integer n, doublereal *x)
{
	// check not two Independents in a row
	assert( ! Recording );

	// make sure the set of independent variables not empty
	assert( n > 0 );

	// check for previously defined function (non-zero pointers)
	if( Xmemory != 0 )
	{	assert( Xmemory != 0 );
		delete [] Xmemory;

		assert( Ymemory != 0 );
		delete [] Ymemory;

		assert( X != 0 );
		delete [] X;

		assert( Y != 0 );
		delete [] Y;
	}

	// allocate memory for independent variable value and differential
	// and store domain value in first column of X
	Xmemory = new double[2 * n];
	X       = new double *[n];
	int j;
	for(j = 0; j < n; j++)
	{	X[j]    = Xmemory + 2 * j;
		X[j][0] = x[j].value();
	}

	// start recording on tape number 1
	Recording = true;
	Domain    = n;
	Range     = 0;
	int tag   = 1;
	int keep  = 1;
	trace_on(tag, keep);

	// declare independent variables
	for(j = 0; j < n; j++)
		x[j] <<= X[j][0];

	return;
}

void Dependent(integer m, doublereal *y)
{
	// check that directly follows call to Independent
	assert( Recording );
	assert( Range == 0 );

	// make sure the set of dependent variables is not empty
	assert( m > 0 );

	// declare the independent variables
	int i;
	double yi;
	for(i = 0; i < m; i++)
		y[i] >>= yi;

	// allocate memory for dependent variable value and differential
	// and store domain value in first column of X
	Ymemory = new double[2 * m];
	Y       = new double *[m];
	for(i = 0; i < m; i++)
		Y[i] = Ymemory + 2 * i;

	// stop the recording
	Range     = m;
	Recording = false;
	trace_off();

}

// linked here from Partial<doublereal>(integer i, integer j)
double Partial(const doublereal &zero, integer i, integer j)
{
	// check that function is defined
	assert( ! Recording );
	assert( Domain > 0 );

	// check dimension of arguments
	int m = Range;
	int n = Domain;
	assert( 0 <= i );
	assert( 0 <= j );
	assert( i < m  );
	assert( j < n );
	assert( zero == doublereal(0) );

	// move data into structure expected by Adolc
	int k;
	for(k = 0; k < Domain; k++)
		X[k][1] = 0.;
	X[j][1] = 1.;

	// compute dy
	int tag  = 1; // tape identifier
	int keep = 1; // keep data in tape
	int d    = 1; // highest derivative degree
	forward(tag, m, n, d, keep, X, Y);

	// return the result
	return Y[i][1];
}

double Value(const doublereal &z)
{	doublereal tmp = z;
	return tmp.value(); 
}

doublereal Sign(const doublereal* x, const doublereal* y)
{	if( *y >= 0. )
		return abs(*x);
	return -abs(*x);
}

} // end f2cad namespace

 // assume that F2CAD_USE_CppAD is defined
# include <f2cad/f2cad.hpp>

namespace {
	
	static CppADvector< CppAD::AD<double> > X;
	static CppAD::ADFun<double> *F = 0;
}

namespace f2cad {

void Independent(integer n, doublereal *x)
{
	// check not two Independents in a row
	assert( X.size() == 0 );

	// make sure the set of independent variables not empty
	assert( n > 0 );

	// check for a previously defined function (nonzero pointer)
	if( F != 0 )
	{	delete F;
		F = 0;
	}

	// move data into structure expected by CppAD
	X.resize(n);
	integer j;
	for(j = 0; j < n; j++)
		X[j] = x[j];

	// declare independent variables
	CppAD::Independent(X);

	// make x equivalent to X
	for(j = 0; j < n; j++)
		x[j] = X[j];

	return;
}

void Dependent(integer m, doublereal *y)
{
	// check that directly follows call to Independent
	assert( X.size() > 0 );
	assert( F == 0 );

	// make sure the set of dependent variables not empty
	assert( m > 0 );

	// move data into structure expected by CppAD
	CppADvector< CppAD::AD<double> > Y(m);
	integer i;
	for(i = 0; i < m; i++)
		Y[i] = y[i];

	// construct the AD function object
	F = new CppAD::ADFun<double>(X, Y);

	// no longer need independent variable information
	X.resize(0);
}


// linked here from Partial<doublereal>(integer i, integer j)
double Partial(const doublereal &zero, integer i, integer j)
{
	// check that directly follows call to Dependent
	assert( X.size() == 0 );
	assert( F != 0 );

	// check dimension of arguments
	integer m = F->Range();
	integer n = F->Domain();
	assert( 0 <= i );
	assert( 0 <= j );
	assert( i < m );
	assert( j < n );
	assert( zero == doublereal(0) );

	// move data into structure expected by CppAD
	CppADvector<double> dx(n);
	integer k;
	for(k = 0; k < n; k++)
		dx[k] = 0.;
	dx[j] = 1.; 

	// Forward mode calculation
	CppADvector<double> dy(m);
	dy = F->Forward(1, dx);

	// return the result
	return dy[i];
}

double Value(const doublereal &z)
{	// convert in a way that works even if we are recording
	return CppAD::Value( Var2Par(z) ); 
}

doublereal Sign(const doublereal* x, const doublereal* y)
{	if( *y >= 0. )
		return abs(*x);
	return -abs(*x);
}

} // end f2cad namespace

 // assume that F2CAD_USE_Fadbad is defined
# include <f2cad/f2cad.hpp>

namespace {
	static bool Recording  = false; // currently recording a function
	static int Domain      = 0;     // domain dimension
	static int Range       = 0;     // range dimension
	static doublereal *Fun = 0;     // currently defined function object
}

namespace f2cad {

void Independent(integer n, doublereal *x)
{
	// check not two Independents in a row
	assert( ! Recording );

	// make sure the set of independent variables not empty
	assert( n > 0 );

	// delete old function object
	if( Fun != 0 )
	{	delete [] Fun;
		Fun = 0;
	}

	// declare independent variables
	int j;
	for(j = 0; j < n; j++)
		x[j].diff(j, n);

	// change to recording state
	Recording = true;
	Domain    = n;
	Range     = 0;

	return;
}

void Dependent(integer m, doublereal *y)
{
	// check that directly follows call to Independent
	assert( Recording );
	assert( Range == 0 );
	assert( Fun == 0 );

	// make sure the set of dependent variables is not empty
	assert( m > 0 );

	// create the function object
	Fun = new doublereal[m];
	int i;
	for(i = 0; i < m; i++)
		Fun[i] = y[i];

	// stop the recording
	Range     = m;
	Recording = false;
}


// linked here from Partial<doublereal>(integer i, integer j)
double Partial(const doublereal &zero, integer i, integer j)
{
	// check that function is defined
	assert( ! Recording );
	assert( Fun != 0 );

	assert( 0 <= i );
	assert( 0 <= j );
	assert( i < Range );
	assert( j < Domain );
	assert( zero == doublereal(0) );

	// compute dy
	return Fun[i].d(j);
}

double Value(const doublereal &z)
{	// value member function is not const
	doublereal tmp = z;
	return tmp.x(); 
}

doublereal Sign(const doublereal* x, const doublereal* y)
{	if( *y >= 0. )
		return abs(*x);
	return -abs(*x);
}

} // end f2cad namespace