Re: Complex Exponential Integral

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

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
> >
> >
> >
> >
> 


-- 
Ethan A Merritt
Department of Biochemistry
University of Washington, Seattle