Re: Complex Exponential Integral

[Ethan A M. <me...@uw...>]

I have just pushed changes to the configure script and to amos_airy,c
so that ZEXINT is used if found during configuration in preference to CEXINT.
I tested configure+build against libraries with or without ZEXINT.

commit dda9e5051d

Since both routines were available for testing, I also sampled the real/imaginary
plane for the difference in output between parallel calls to CEXINT and ZEXINT
for orders 2,3,4, and 5.
The maximum difference found was on the order of    3 x 10^{-15} 
Pure noise level for double precision would be about  2 x 10^{-16}
Pure noise level for single precision would be about    1 x 10^{-7}.

cheers,
	Ethan

On Wednesday, 9 July 2025 23:55:03 PDT Ethan A Merritt wrote:
> On Wednesday, 9 July 2025 22:12:34 PDT Erik Luijten wrote:
> > I think we are mostly in agreement; I would be happy with "configure" being
> > able to detect ZEXINT and then use the call that I implemented in
> > amos_airy.c (as attached to my earlier email).
> 
> Right.  I'll work on a modification to the configure script.
> 
> > Just for completeness: My concerns with the current situation are that (1)
> > it is not documented (unless I overlooked it) that CEXINT must be compiled
> > in a way that upgrades all floats to double (plus accompanying changes to
> > i1mach/d1mach); 
> 
> [shrug]  I think it's outside the scope of gnuplot documentation to instruct
> people on how to build a math library.   I do bundle build instructions and a
> Makefile together with the libamos source files if anyone asks me for a copy.
> I put this all together quite a while ago, but at the time I found I needed to
> tell the compiler to make REAL*8 and DOUBLE*8 the default for all the
> Amos routines, not just CEXINT.  Makefile attached. 
> 
> > (2) *generally* there is no guarantee that such an
> > "upgrade" would yield the desired precision. I acknowledge that for CEXINT
> > it may work fine, so this is primarily an aesthetic objection.
> 
> I have just done a line-by-line comparison of the ZEXINT and CEXINT code.
> Except for the D vs E declaration of Fortran constant values which I noted
> before, the only differences are the packing/unpacking of real and imaginary
> components into a single COMPLEX variable vs two parallel DOUBLE variables.
> So yes, they two versions of the algorithm have the same precision when
> compiled with the same compiler flags. 
> 
> Note: The 1990 Amos paper documents the accuracy of the algorithm as
> "18 digits" but notes that this may not be achieved on the hardware of
> "some machines".  And indeed this is higher precision than guaranteed
> by IEEE 754 and 64-bits. I suspect anyone who cares about the difference
> between 16-bit and 18-bit accuracy for expint() is not relying on gnuplot
> to provide it :-)
> 
> Digression:   When I modified gnuplot to use 64-bit integers I also looked
> into using gcc extensions so that the floating point calculations could use
> 80-bits rather than 64-bits on architectures that support it.  That turned
> out to be a morass; I dropped the idea pretty fast. 
> 
> > Please let me know once you are ready to let me test the new version of
> > configure.
> 
> Sure.  I'll probably break out and add ZEXINT itself to my library collection
> first so I will be able to make sure the whole thing works here.   
> 
>   cheers
> 	Ethan
> 
> > Best,
> > 
> > Erik
> > 
> > 
> > On Wed, Jul 9, 2025 at 9:30 PM Ethan A Merritt <me...@uw...> wrote:
> > 
> > > On Wednesday, 9 July 2025 06:25:07 PDT Erik Luijten wrote:
> > > >
> > > > - That being said, I think it is not ideal that this configuration option
> > > > then uses a version of cexint that must be modified by the user. Instead,
> > > > the version of amos_airy.c that I sent yesterday uses an unmodified
> > > ZEXINT
> > > > as distributed in Netlib, minimizing confusion (and effort) for anyone
> > > who
> > > > wishes to incorporate it when they compile gnuplot. You clearly
> > > understand
> > > > the potential issues of upgrading a single-precision function to double
> > > > precision, but not everyone does and it also would require the user to
> > > > address a matter that was resolved by the original author from the
> > > outset.
> > >
> > > I am not following your thought there.  What matter is that?
> > >
> > > The biggest change/cleanup/upgrade/whatever is that the original code
> > > was written in an era when hardware floating support varied widely.
> > > The r1mach/d1mach per-installation customization had to deal with
> > > all sorts of setups, whereas now hardware support for IEEE 754 floating
> > > point is nearly universal.
> > >
> > > As described in the 1990 paper, the separate double-precision routine was
> > > needed because at that time Fortran did not support double precision
> > > complex
> > > as an intrinsic type. Since then (Fortran90? Fortran95?μ) this limitation
> > > went away, and so long as the algorithm itself provides sufficient
> > > precision
> > > the use of single vs double precision can be made a compile-time option.
> > >
> > > I do not recall having to change anything in the code to CEXINT itself.
> > > I only replaced the separate machine-dependent routines in favor of
> > > letting the compiler pick up the relevant constants from the system
> > > environment it was compiled on, and the compiler flags used to build it.
> > >
> > > > - By the way, you mention that you found an older version of CEXINT, but
> > > > the text you quote ("REVISION DATE  870515", etc.) is identical in the
> > > ACM
> > > > version https://urldefense.com/v3/__https://www.netlib.org/toms-2014-06-10/683__;!!K-Hz7m0Vt54!iwiffB-qjYH_1hni74sIr0ktLbCn0P2jQlKbCMP8UgcD3GZ7yPj0ouUn2UsTg7d08Hwtnmzu50gBQtiWqElv$  ; there are just
> > > other
> > > > subroutines preceding it in the same file. In fact, I do not think that
> > > the
> > > > copyright situation favors CEXINT; the two functions were created at the
> > > > same time under the same circumstances (ZEXINT also carries the 1987
> > > date).
> > >
> > > Sure.  Both routines were written at the same time.  When the source code
> > > was placed in the TOMS archive someone added the comment lines at the top
> > > citing the 1990 paper.  My only point was that copies of the source code
> > > existed
> > > that predated the TOMS publication thus avoided any implied claim of
> > > license
> > > requirements imposed by the journal.
> > >
> > > > - Regarding concrete compilation: I compile ZEXINT and the other
> > > functions
> > > > it uses into a separate library (unlike you, I use i1mach and d1mach, but
> > > > rely on those being provided in libopenspecfun). I found that if I
> > > specify
> > > > this library in LIBAMOS_LDFLAGS, configure will detect it. Practically
> > > > speaking, I make this a static library so that users of the macOS binary
> > > do
> > > > not need to worry about it (I distribute gnuplot as a single static
> > > > executable).
> > > >
> > > > Depending on your view, the only question is then whether we switch from
> > > > CEXINT to ZEXINT (which would also require a small change in
> > > "configure").
> > >
> > > I wouldn't "switch".  I would have configure search for both and use
> > > whichever
> > > it found first.  I don't care much which one it finds first, unless as you
> > > say there
> > > is the possibility of conflicting symbols in another library.  If I
> > > understand your
> > > plan correctly, that would not be a problem for either of us.
> > >
> > > The autoconf scripts are not my strong point (to say the least), but I'll
> > > see
> > > what I can come up with.
> > >
> > > > Apologies for reviving this old topic!
> > >
> > > No problem - a trip down memory lane.
> > >
> > > By the way, the Julia language project uses a reimplemention in c99 for
> > > the Amos library functions including expint().   I did not pursue using it
> > > because unlike openspecfun it would not be available to users as a standard
> > > linux distribution package.
> > > But if you're already building static libraries yourself, and not under
> > > linux,
> > > that might be of interest to you as an alternative.
> > >
> > >   cheers,
> > >         Ethan
> > >
> > > >
> > > > Erik
> > > >
> > > > On Tue, Jul 8, 2025 at 11:52 PM Ethan A Merritt <me...@uw...> wrote:
> > > >
> > > > > On Tuesday, 8 July 2025 18:55:03 PDT Erik Luijten wrote:
> > > > > > Hi,
> > > > > >
> > > > > > I found the code on
> > > https://urldefense.com/v3/__https://dl.acm.org/doi/abs/10.1145/78928.78934__;!!K-Hz7m0Vt54!jVJeG1oeF8DBF4-meLyPbDTC2nBP2lu0ChoIkfgKeOZGsqFA4LqPWSqtJ5ZK0ZYv1jwibOIAGDbUEM2Nmkq6$
> > > ; the
> > > > > > supplementary material contains all the subroutines concatenated
> > > into one
> > > > > > file.
> > > > > > It is also available on
> > > https://urldefense.com/v3/__https://www.netlib.org/toms-2014-06-10/__;!!K-Hz7m0Vt54!jVJeG1oeF8DBF4-meLyPbDTC2nBP2lu0ChoIkfgKeOZGsqFA4LqPWSqtJ5ZK0ZYv1jwibOIAGDbUENvigquW$
> > > > > (scroll to
> > > > > > algorithm 683).
> > > > > > The fact that it is on netlib tells me it is likely free to use.
> > > > >
> > > > > At the top of that page is a notice:
> > > > >
> > > > >     Use of ACM Algorithms is subject to the ACM Software Copyright and
> > > > >     License Agreement
> > > > >
> > > > > > However, the ACM TOMS license agreement is here:
> > > > > >
> > > https://urldefense.com/v3/__https://www.acm.org/publications/policies/software-copyright-notice__;!!K-Hz7m0Vt54!jVJeG1oeF8DBF4-meLyPbDTC2nBP2lu0ChoIkfgKeOZGsqFA4LqPWSqtJ5ZK0ZYv1jwibOIAGDbUEHrkADU9$
> > > .
> > > > > > I am not sure if this conflicts with Gnuplot's license?
> > > > >
> > > > > The issue is not that it conflicts with gnuplot's license.
> > > > > The issue is that Debian, for example, was concerned that inclusion
> > > > > with gnuplot would make the whole package non-free.
> > > > >
> > > > > Thanks for the link to the ACM site.  They have very much changed
> > > > > their policy since I first looked into this.  And it seems that for
> > > > > software
> > > > > published since 2013 the ACM does not claim any rights.  That's good.
> > > > > Unfortunately it doesn't help.  For older software (the Amos routines
> > > > > are from the 1980s) they state that the terms of the
> > > > > earlier ACM license still apply, which restricts commercial use.
> > > > > If that were true it would make the whole set of Amos routines
> > > > > [as obtained from them] non-free, but see below.
> > > > >
> > > > > > Whether a more permissive license applies to cexint_.f  than to
> > > > > > zexint_.f is also unclear to me.
> > > > >
> > > > > It is my belief that ACM was simplly wrong in ever claiming the right
> > > > > to license or restrict any of these routines.  The were produced at
> > > > > Sandia National Laboratories under US government contract, and as
> > > > > I understand it (I am not a lawyer and all that) that makes them
> > > > > "Works Created by the United States Government" and as such not
> > > > > subject to copyright.
> > > > >
> > > > > But that's only me.  I don't have any sway over what Debian or anyone
> > > > > else considers to be "non-free".  Based on past conversations my
> > > > > understanding is that Debian considers the TOMS license to poison
> > > > > any project that includes it as non-free.  Those conversations were
> > > > > many years ago, so I suppose the issue could be revisited.
> > > > >
> > > > > > An awkward point, to the best of my understanding, is that an
> > > unmodified
> > > > > > version of cexint will not work with the Gnuplot source code, as it
> > > > > expects
> > > > > > single-precision arguments. Yet, changing cexint is not
> > > > > > completely straightforward: Comparing cexint and zexint (and the
> > > > > > accompanying routines cexent and zexent), I noticed the more subtle
> > > point
> > > > > > that cexint uses r1mach.f and zexint uses d1mach.f; moreover, they
> > > use
> > > > > > different constants from i1mach to establish precision (elements
> > > 12,13
> > > > > vs.
> > > > > > 15,16). Also, in zexint all constants are replaced by their
> > > > > > double-precision counterparts.
> > > > >
> > > > > This gets complicated.
> > > > >
> > > > > I have downloaded and had a quick look at the code available through
> > > the
> > > > > netlib site.  It is apparently newer than the code I have here. It
> > > cites
> > > > > the 1990
> > > > > TOMS paper at the top and as you say it concatenates many smaller
> > > routines,
> > > > > whereas the code I have claims to be from 1987 and starts
> > > > > C***BEGIN PROLOGUE  CEXINT
> > > > > C***DATE WRITTEN   870515   (YYMMDD)
> > > > > C***REVISION DATE  870515   (YYMMDD)
> > > > > C***CATEGORY NO.  B5E
> > > > > C***KEYWORDS  EXPONENTIAL INTEGRALS, SINE INTEGRAL, COSINE INTEGRAL
> > > > > C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
> > > > > C***PURPOSE  TO COMPUTE EXPONENTIAL INTEGRALS OF A COMPLEX ARGUMENT
> > > > >
> > > > > There are some code differences elsewhere, but I did not look in detail
> > > > > and don't have anything useful to say about them right now.
> > > > >
> > > > > Beyond that, to the extent that there are precision-sensitive
> > > constants in
> > > > > the
> > > > > code the numerical values seem identical, but they are declared, for
> > > > > example
> > > > > -5.77215664901532861E-01   vs   -5.77215664901532861D-01
> > > > > As I recall in 1987 that difference was necessary to distinguish
> > > single-
> > > > > and
> > > > > double- precision constants, but I don't think current Fortran pays any
> > > > > attention to it when using the constant in a double-precision context.
> > > > > I compile with
> > > > >      gfortran -fdefault-real-8 -fdefault-double-8
> > > > > I.e. everything should be double precision by default.
> > > > > So although you describe cexint() as a single-precision routine,
> > > > > that isn't true when I build it here.  Of course it's always possible
> > > > > I missed something in the code that loses precision, but given that
> > > > > it iterates until convergence I think that would not even matter.
> > > > >
> > > > > > For all these reasons, I would argue that if copyright restrictions
> > > are
> > > > > the
> > > > > > same for cexint and zexint, it would be better to use the latter - a
> > > user
> > > > > > could then directly download this subroutine and use it without
> > > > > > modifications.
> > > > >
> > > > > You may have a point about an end-user who builds from source.
> > > > > But you've been building and making available binaries for MacOS,
> > > right?
> > > > > I think that moves any legal question to "distribution" rather than
> > > "use".
> > > > > You may well agree with me that the code you want to include is not
> > > subject
> > > > > to the TOMS license even though TOMS seems to say it is.  But you
> > > should
> > > > > at least be aware that the question has been raised in the past.
> > > > >
> > > > > > Lastly, a tricky point that is easily overlooked is that i1mach,
> > > r1mach,
> > > > > > d1mach need to be configured for x86. Moreover, since openspecfun is
> > > > > built
> > > > > > as a library, it already contains i1mach and d1mach, opening the
> > > door for
> > > > > > conflicting functions if one gathers cexint/zexint into a separate
> > > > > > library...
> > > > >
> > > > > As to r1mach.f and d1mach.f,  I am not using either of them here.
> > > > > I replaced these with a C-language data block (attached) that I use
> > > > > to build libamos from the Sandia sources for all the Amos routines.
> > > > > I link against that rather than against libopenspecfun.
> > > > >
> > > > > I sent my collection of libamos source to Tatsuro Matsuoka in case
> > > > > he wanted to use them for building Windows packages.
> > > > > I even  considered adding source for all the libamos routines as a
> > > > > separate subdirectory in the gnuplot distribution package.
> > > > > But then I found that many/most linux distros  already provide
> > > > > libopenspecfun so I decided against complicating the gnuplot packaging
> > > > > by having it supply a whole separate library for the sake of a single
> > > > > routine.
> > > > >
> > > > > > I am attaching my modified version of amos_airy.c, in case you are
> > > swayed
> > > > > > by my argument...
> > > > >
> > > > > Are you thinking that the user would build zexint.o as a separate
> > > object
> > > > > module and then add it to the Makefile somehow?
> > > > >
> > > > > I would prefer to make it a configuration option. Not sure how best
> > > > > to capture a single extra source file whose provenance and dependencies
> > > > > we don't control.  I think it would be easier to require that it has
> > > > > already
> > > > > been wrapped into a shared library.
> > > > > Or hey, another alternative is to provide a template for making it
> > > > > a plugin.  The distribution package currently has a demo plugin for
> > > > > uigamma.  Providing a plugin for a Fortran routine might be of
> > > > > interest even apart from filling a gap in special function coverage.
> > > > >
> > > > >       Ethan
> > > > >
> > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > On Tue, Jul 8, 2025 at 7:26 PM Ethan A Merritt <me...@uw...>
> > > wrote:
> > > > > >
> > > > > > > On Tuesday, 8 July 2025 14:21:01 PDT Erik Luijten wrote:
> > > > > > > > Hi Ethan,
> > > > > > > >
> > > > > > > > In amos_airy.c, you write:
> > > > > > > >
> > > > > > > >  * The cexint_() routine is also from the AMOS collection, but
> > > is not
> > > > > > > part
> > > > > > > >  * of the Bessel function subset and not included in
> > > libopenspecfun.
> > > > > > > >  * I could not find source for a version of the source that
> > > splits
> > > > > the
> > > > > > > >  * real and imaginary parts of the arguments as with the zxxxx.f
> > > > > Bessel
> > > > > > > >  * routines; this one accepts and returns a CMPLX argument rather
> > > > > than
> > > > > > > > separate
> > > > > > > >  * real and imaginary parts.
> > > > > > > >  *             D. E. Amos, ALGORITHM 683  A Portable FORTRAN
> > > > > Subroutine
> > > > > > > >  *             for Exponential Integrals of a Complex Argument
> > > > > > > >  *             ACM Trans. Math. Software 16:178-182 (1990)
> > > > > > > >
> > > > > > > > This confuses me for two reasons:
> > > > > > > > - You invoke CEXINT as a double precision function, even though
> > > this
> > > > > > > > Fortran function is single precision.
> > > > > > >
> > > > > > > > - The algorithm 683 you quote DOES provide a separate double
> > > > > precision
> > > > > > > > function, ZEXINT. Moreover, that function does exactly what you
> > > were
> > > > > > > > looking for, namely taking and returning separate real and
> > > imaginary
> > > > > > > parts.
> > > > > > >
> > > > > > > There is some history behind this due to IMHO over-zealous, or at
> > > least
> > > > > > > overly pessimistic, concern about license contamination with
> > > non-free
> > > > > code.
> > > > > > > I cited a publication in ACM/TMS that describes the implementation.
> > > > > > > However, people have assserted that code obtained through the ACM
> > > > > > > journal is non-free because the journal claims copyright.
> > > > > > > I believe this to be factually incorrect; the published article
> > > may be
> > > > > > > subject
> > > > > > > to copyright but not the code it describes. In any case, rather
> > > than
> > > > > pursue
> > > > > > > the argument I instead used copies of the code obtained from US
> > > > > national
> > > > > > > lab repositories.  I found source there for cexint.f but not for
> > > > > zexint.f,
> > > > > > > so that is what I use myself to build libamos and what I invoked
> > > from
> > > > > the
> > > > > > > gnuplot code.
> > > > > > >
> > > > > > > Even so that was apparently not sufficient by itself to convince
> > > > > hard-core
> > > > > > > "free licenses only" gatekeepers.
> > > > > > > There is an explicit statement of the [non]copyright status of the
> > > code
> > > > > > > in the cited 1985 Sandia National Laboratory report, but that
> > > report
> > > > > > > predated development of the expint routines (1990?) and I did not
> > > find
> > > > > a
> > > > > > > separate Sandia report that mentions them specifically.
> > > > > > > That is probably why neither routine is included in libopenspec,
> > > > > > > and why inclusion in gnuplot is a configuration option.
> > > > > > >
> > > > > > > > I changed f_amos_cexint to use ZEXINT (just for my own
> > > edification)
> > > > > and
> > > > > > > it
> > > > > > > > works properly. However, before submitting it, I wonder what am I
> > > > > > > missing?
> > > > > > >
> > > > > > > Where did you obtain a copy of the source?
> > > > > > > Can you document the relevant claim or non-claim of copyright or
> > > > > > > licensing?
> > > > > > >
> > > > > > > And, out of curiosity, is there anything in the zexint.f code that
> > > > > would
> > > > > > > differentiate it from a single-precision version other than
> > > declaration
> > > > > > > of the variables as Fortran DOUBLE PRECISION?
> > > > > > >
> > > > > > >         Ethan
> > > > > > >
> > > > > > >
> > > > > > > >
> > > > > > > > Thank you,
> > > > > > > >
> > > > > > > > Erik
> > >
> > >
> > >
> > >
> > 
> 
> 
>