Re: Complex Exponential Integral

This is a multi-part message in MIME format.

--nextPart3740467.HovnAMPojK
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain; charset="UTF-8"

On Tuesday, 8 July 2025 18:55:03 PDT Erik Luijten wrote:
> Hi,
>=20
> I found the code on https://dl.acm.org/doi/abs/10.1145/78928.78934 ; the
> supplementary material contains all the subroutines concatenated into one
> file.
> It is also available on https://www.netlib.org/toms-2014-06-10/  (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://www.acm.org/publications/policies/software-copyright-notice.
> I am not sure if this conflicts with Gnuplot's license?=20

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.  =20

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 expec=
ts
> 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 th=
e 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 exampl=
e=20
=2D5.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=20
     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 t=
he
> 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.=20
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 routi=
ne.

> 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


>=20
>=20
>=20
> On Tue, Jul 8, 2025 at 7:26=E2=80=AFPM Ethan A Merritt <[email protected]> w=
rote:
>=20
> > 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 Bess=
el
> > >  * routines; this one accepts and returns a CMPLX argument rather than
> > > separate
> > >  * real and imaginary parts.
> > >  *             D. E. Amos, ALGORITHM 683  A Portable FORTRAN Subrouti=
ne
> > >  *             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 c=
ode.
> > 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 pu=
rsue
> > 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=
=2Ef,
> > so that is what I use myself to build libamos and what I invoked from t=
he
> > gnuplot code.
> >
> > Even so that was apparently not sufficient by itself to convince hard-c=
ore
> > "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) a=
nd
> > 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
> > >
> >
> >
> >
> >
>=20


--nextPart3740467.HovnAMPojK
Content-Disposition: attachment; filename="d1mach.c"
Content-Transfer-Encoding: 7Bit
Content-Type: text/x-csrc; charset="x-UTF_8J"; name="d1mach.c"

/* Ethan A Merritt - Aug 2012
 * This version of D1MACH is for use with
 * gfortran -fdefault-real-8 -fdefault-double-8
 * i.e. 8 byte REAL and DOUBLE but 4 byte integers
 */
#include <stdint.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
double d1mach_(int32_t *i)
{
	switch(*i){
	  case 1: return DBL_MIN;
	  case 2: return DBL_MAX;
	  case 3: return DBL_EPSILON/FLT_RADIX;
	  case 4: return DBL_EPSILON;
	  case 5: return log10((double)FLT_RADIX);
	  }
	fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
}

--nextPart3740467.HovnAMPojK
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline


--nextPart3740467.HovnAMPojK
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
gnuplot-beta mailing list
[email protected]
Membership management via: https://lists.sourceforge.net/lists/listinfo/gnuplot-beta

--nextPart3740467.HovnAMPojK--