hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

convert docs to sphinx

Browse source
This commit is contained in:
Fredrik Johansson 2012-10-11 13:09:44 +02:00
parent ef60ec1779
commit 90660611af
21 changed files with 2275 additions and 2079 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
README
View file

@ -4,7 +4,7 @@ Planned features include complex numbers and fast computation of special
functions. The project might get merged back into FLINT when the code
stabilizes in the future.
For documentation and additional details, see: http://fredrikj.net/arb/doc.html
For documentation and additional details, see: http://fredrikj.net/arb/
Development updates are occasionally posted to: http://fredrikj.net/blog/

153
doc/Makefile Normal file
View file

@ -0,0 +1,153 @@
# Makefile for Sphinx documentation
#
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
PAPER =
BUILDDIR = build
# Internal variables.
PAPEROPT_a4 = -D latex_paper_size=a4
PAPEROPT_letter = -D latex_paper_size=letter
ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
# the i18n builder cannot share the environment and doctrees with the others
I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext
help:
@echo "Please use \`make <target>' where <target> is one of"
@echo " html to make standalone HTML files"
@echo " dirhtml to make HTML files named index.html in directories"
@echo " singlehtml to make a single large HTML file"
@echo " pickle to make pickle files"
@echo " json to make JSON files"
@echo " htmlhelp to make HTML files and a HTML help project"
@echo " qthelp to make HTML files and a qthelp project"
@echo " devhelp to make HTML files and a Devhelp project"
@echo " epub to make an epub"
@echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
@echo " latexpdf to make LaTeX files and run them through pdflatex"
@echo " text to make text files"
@echo " man to make manual pages"
@echo " texinfo to make Texinfo files"
@echo " info to make Texinfo files and run them through makeinfo"
@echo " gettext to make PO message catalogs"
@echo " changes to make an overview of all changed/added/deprecated items"
@echo " linkcheck to check all external links for integrity"
@echo " doctest to run all doctests embedded in the documentation (if enabled)"
clean:
-rm -rf $(BUILDDIR)/*
html:
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
dirhtml:
$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
singlehtml:
$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
@echo
@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
pickle:
$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
@echo
@echo "Build finished; now you can process the pickle files."
json:
$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
@echo
@echo "Build finished; now you can process the JSON files."
htmlhelp:
$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
@echo
@echo "Build finished; now you can run HTML Help Workshop with the" \
".hhp project file in $(BUILDDIR)/htmlhelp."
qthelp:
$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
@echo
@echo "Build finished; now you can run "qcollectiongenerator" with the" \
".qhcp project file in $(BUILDDIR)/qthelp, like this:"
@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/Arb.qhcp"
@echo "To view the help file:"
@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/Arb.qhc"
devhelp:
$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
@echo
@echo "Build finished."
@echo "To view the help file:"
@echo "# mkdir -p $$HOME/.local/share/devhelp/Arb"
@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/Arb"
@echo "# devhelp"
epub:
$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
@echo
@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
latex:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo
@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@echo "Run \`make' in that directory to run these through (pdf)latex" \
"(use \`make latexpdf' here to do that automatically)."
latexpdf:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo "Running LaTeX files through pdflatex..."
$(MAKE) -C $(BUILDDIR)/latex all-pdf
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
text:
$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
@echo
@echo "Build finished. The text files are in $(BUILDDIR)/text."
man:
$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
@echo
@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
texinfo:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo
@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
@echo "Run \`make' in that directory to run these through makeinfo" \
"(use \`make info' here to do that automatically)."
info:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo "Running Texinfo files through makeinfo..."
make -C $(BUILDDIR)/texinfo info
@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
gettext:
$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
@echo
@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
changes:
$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
@echo
@echo "The overview file is in $(BUILDDIR)/changes."
linkcheck:
$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
@echo
@echo "Link check complete; look for any errors in the above output " \
"or in $(BUILDDIR)/linkcheck/output.txt."
doctest:
$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
@echo "Testing of doctests in the sources finished, look at the " \
"results in $(BUILDDIR)/doctest/output.txt."

8
doc/credits.txt
View file

@ -1,8 +0,0 @@
<p>Arb is licensed GNU General Public License version 2, or any later version.</p>
<p>Fredrik's work on Arb is supported by Austrian Science Fund FWF Grant Y464-N18
(Fast Computer Algebra for Special Functions).</p>
<p>Arb includes code written by Bill Hart and Sebastian Pancratz
for FLINT (also licensed GPL 2.0+).</p>

181
doc/doc.py
View file

@ -1,181 +0,0 @@
import datetime
import time
tz = time.tzname[0]
now = datetime.datetime.now()
intro = r"""
<html>
<head>
<title>Arb documentation</title>
<style>
body { font-family: sans-serif; font-size: 90%; }
dt { color:#348; font-weight: bold; }
dd p { margin-top: 0.4em; }
</style>
<script type='text/x-mathjax-config'>
MathJax.Hub.Config({
extensions: ["tex2jax.js"],
jax: ["input/TeX", "output/HTML-CSS"],
MMLorHTML: { prefer: "HTML" },
tex2jax: {
inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
processEscapes: true
},
"HTML-CSS": { availableFonts: ["TeX"] }
});</script>
<script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
<script type='text/javascript'>
newContainer = document.createElement('span');
newContainer.style.setProperty("display","none","");
newNode = document.createElement('script');
newNode.type = "math/tex";
newNode.innerHTML = '\\newcommand{\\NN}{\\mathbb{N}}\n\\newcommand{\\ZZ}{\\mathbb{Z}}\n\\newcommand{\\QQ}{\\mathbb{Q}}\n\\newcommand{\\FF}{\\mathbb{F}}\n\\newcommand{\\KK}{\\mathbb{K}}\n\\newcommand{\\RR}{\\mathbb{R}}\n\\newcommand{\\CC}{\\mathbb{C}}\n\\newcommand{\\abs}[1]{\\left|#1\\right|}';
newContainer.appendChild(newNode);
document.body.insertBefore(newContainer,document.body.firstChild);
</script>
</head>
<body>
<h1>Arb documentation</h1>
<p><i>Last updated: %NOW%</i></p>
<h2>Contents</h2>
"""
intro = intro.replace("%NOW%", "%i-%02i-%02i %02i:%02i:%02i %s" % (now.year, now.month, now.day, now.hour, now.minute, now.second, tz))
outro = r"""
<p>Wishes or bug reports? Send any comments to <a href="mailto:fredrik.johansson@gmail.com">fredrik.johansson@gmail.com</a></p>
</body>
</html>
"""
def section_key(title, section):
def stringify(s):
t = ""
for c in s:
if c.isalpha():
t += c
else:
t += "-"
return t
return stringify(title) + "-" + stringify(section)
def parse(string, markup):
if not markup:
return string
lines = string.splitlines()
# Collect paragraphs
chunks = []
chunk = []
for line in lines:
if not line.strip():
if chunk:
chunks.append(chunk)
chunk = []
else:
chunk.append(line)
doc = []
heading = ""
description = []
in_heading = False
in_description = False
chunks.append(["***"])
for chunk in chunks:
if len(chunk) == 1 and chunk[0].startswith("***"):
if in_heading:
in_heading = False
doc.append(("H", heading.lstrip()))
heading = ""
else:
in_heading = True
if in_description:
doc.append(("DESCR", description))
description = []
in_description = False
else:
if chunk[0].startswith(" "):
if in_heading:
heading += " ".join(chunk)
else:
in_description = True
description.append(" ".join(chunk))
else:
if in_description:
doc.append(("DESCR", description))
description = []
in_description = False
doc.append(("DEF", " ".join(chunk)))
return doc
def write(docs, file):
file.write(intro)
file.write("<ul>\n");
# Write TOC
for doc, title in docs:
file.write("<li>\n");
file.write('<a href="#%s">%s</a>\n' % (section_key(title, ""), title));
file.write("<ul>\n");
if isinstance(doc, list):
for entry in doc:
if entry[0] == "H":
file.write('<li><a href="#%s">%s</a></li>\n' % (section_key(title, entry[1]), entry[1]))
file.write("</ul></li>\n");
file.write("</ul>\n");
file.write("<hr />\n")
# Write each section
for doc, title in docs:
in_list = False
file.write('<h2><a name="%s">%s</a></h2>\n' % (section_key(title, ""), title));
if isinstance(doc, list):
for entry in doc:
if entry[0] == "H":
if in_list:
in_list = False
file.write("</dl>\n")
file.write('<a name="%s"><h3>%s</h3></a>\n' % (section_key(title, entry[1]), entry[1]))
if entry[0] == "DEF":
if not in_list:
in_list = True
file.write("<dl>\n")
#file.write("<dt><tt>%s</tt></dt>\n" % entry[1])
file.write("<dt>%s</dt>\n" % entry[1])
if entry[0] == "DESCR":
if in_list:
file.write("<dd>\n")
for para in entry[1]:
file.write("<p>%s</p>\n" % para)
if in_list:
file.write("</dd>\n")
if in_list:
file.write("</dl>\n")
else:
file.write(doc)
file.write("<hr />\n")
file.write(outro)
docs = [
("introduction.txt", "Introduction", False),
("setup.txt", "Setup", False),
("fmpr.txt", "fmpr.h (floating-point arithmetic)", True),
("fmprb.txt", "fmprb.h (real ball arithmetic)", True),
("fmprb_poly.txt", "fmprb_poly.h (polynomials of real balls)", True),
("fmprb_mat.txt", "fmprb_mat.h (matrices of real balls)", True),
("history.txt", "History", False),
("credits.txt", "Credits", False),
]
write([(parse(open(doc).read(), markup), title) for \
(doc, title, markup) in docs], open("doc.html", "w"))

491
doc/fmpr.txt
View file

@ -1,491 +0,0 @@
A variable of type <tt>fmpr_t</tt> holds an arbitrary-precision binary
floating-point number, i.e. a rational number of the form
$x \times 2^y$ where $x, y \in \mathbb{Z}$ and $x$ is odd;
or one of the special values zero, plus infinity, minus infinity,
or NaN (not-a-number).
The component $x$ is called the <i>mantissa</i>, and $y$ is called the
<i>exponent</i>. Note that this is just one among many possible
conventions: the mantissa (alternatively <i>significand</i>) is
sometimes viewed as a fraction in the interval $[1/2, 1)$, with the
exponent pointing to the position above the top bit rather than the
position of the bottom bit, and with a separate sign.
The conventions for special values largely follow those of the
IEEE floating-point standard. At the moment, there is no support
for negative zero, unsigned infinity, or a NaN with a payload, though
some these might be added in the future.
An <tt>fmpr</tt> number is exact and has no inherent "accuracy". We
use the term <i>precision</i> to denote either the target precision of
an operation, or the bit size of a mantissa (which in general is
unrelated to the "accuracy" of the number: for example, the
floating-point value 1 has a precision of 1 bit in this sense and is
simultaneously an infinitely accurate approximation of the
integer 1 and a 2-bit accurate approximation of
$\sqrt 2 = 1.011010100\ldots_2$).
Except where otherwise noted, the output of an operation is the
floating-point number obtained by taking the inputs as exact numbers,
in principle carrying out the operation exactly, and rounding the
resulting real number to the nearest representable floating-point
number whose mantissa has at most the specified number of bits, in
the specified direction of rounding. Some operations are always
or optionally done exactly.
*******************************************************************************
Types, macros and constants
*******************************************************************************
fmpr_struct
An <tt>fmpr_struct</tt> holds a mantissa and an exponent.
If the mantissa and exponent are sufficiently small, their values are
stored as immediate values in the <tt>fmpr_struct</tt>; large values are
represented by pointers to heap-allocated arbitrary-precision integers.
Currently, both the mantissa and exponent are implemented using
the FLINT <tt>fmpz</tt> type. Special values are currently encoded
by the mantissa being set to zero.
fmpr_t
An <tt>fmpr_t</tt> is defined as an array of length one of type
<tt>fmpr_struct</tt>, permitting an <tt>fmpr_t</tt> to be passed by
reference.
fmpr_rnd_t
Specifies the rounding mode for the result of an approximate operation.
FMPR_RND_NEAREST
Specifies that the result of an operation should be rounded to the
nearest representable number, rounding to an odd mantissa if there is a tie
between two values. Note: the code for this rounding mode is currently
not implemented.
FMPR_RND_DOWN
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards zero.
FMPR_RND_UP
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction away from zero.
FMPR_RND_FLOOR
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards minus infinity.
FMPR_RND_CEIL
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards plus infinity.
FMPR_PREC_EXACT
If passed as the precision parameter to a function, indicates that no
rounding is to be performed. This must only be used when it is known
that the result of the operation can be represented exactly and fits
in memory (the typical use case is working with values small integers).
Note that, for example, adding two numbers whose exponents are far
apart can easily produce an exact result that is far too large to
store in memory.
*******************************************************************************
Memory management
*******************************************************************************
void fmpr_init(fmpr_t x)
Initializes the variable x for use. Its value is set to zero.
void fmpr_clear(fmpr_t x)
Clears the variable x, freeing or recycling its allocated memory.
*******************************************************************************
Special values
*******************************************************************************
void fmpr_zero(fmpr_t x)
void fmpr_one(fmpr_t x)
void fmpr_pos_inf(fmpr_t x)
void fmpr_neg_inf(fmpr_t x)
void fmpr_nan(fmpr_t x)
Sets x respectively to 0, 1, $+\infty$, $-\infty$, NaN.
int fmpr_is_zero(const fmpr_t x)
int fmpr_is_one(const fmpr_t x)
int fmpr_is_pos_inf(const fmpr_t x)
int fmpr_is_neg_inf(const fmpr_t x)
int fmpr_is_nan(const fmpr_t x)
Returns nonzero iff x respectively equals
0, 1, $+\infty$, $-\infty$, NaN.
int fmpr_is_inf(const fmpr_t x)
Returns nonzero iff x equals either $+\infty$ or $-\infty$.
int fmpr_is_normal(const fmpr_t x)
Returns nonzero iff x is a finite, nonzero floating-point value, i.e.
not one of the special values 0, $+\infty$, $-\infty$, NaN.
int fmpr_is_special(const fmpr_t x)
Returns nonzero iff x is one of the special values
0, $+\infty$, $-\infty$, NaN, i.e. not a finite, nonzero
floating-point value.
*******************************************************************************
Assignment and rounding
*******************************************************************************
long _fmpr_normalise(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd)
Rounds the mantissa and exponent in-place.
void fmpr_set(fmpr_t y, const fmpr_t x)
Sets y to a copy of x.
long fmpr_set_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets y to a copy of x rounded in the direction specified by rnd to the
number of bits specified by prec.
void fmpr_set_error_result(fmpr_t err, const fmpr_t result, long rret)
Given the return value rret and output variable result from a
function performing a rounding (e.g. fmpr_set_round or fmpr_add), sets
err to a bound for the absolute error.
void fmpr_add_error_result(fmpr_t err, const fmpr_t err_in,
const fmpr_t result, long rret, long prec, fmpr_rnd_t rnd)
Like fmpr_set_error_result, but adds err_in to the error.
*******************************************************************************
Comparisons
*******************************************************************************
int fmpr_equal(const fmpr_t x, const fmpr_t y)
Returns nonzero iff x and y are exactly equal. This function does
not treat NaN specially, i.e. NaN compares as equal to itself.
int fmpr_cmp(const fmpr_t x, const fmpr_t y)
Returns negative, zero, or positive, depending on whether x is respectively
smaller, equal, or greater compared to y. Comparison with NaN is undefined.
int fmpr_cmpabs(const fmpr_t x, const fmpr_t y)
Compares the absolute values of x and y.
int fmpr_sgn(const fmpr_t x)
Returns $-1$, $0$ or $+1$ according to the sign of x. The sign
of NaN is undefined.
*******************************************************************************
Random number generation
*******************************************************************************
void fmpr_randtest(fmpr_t x, flint_rand_t state, long bits, long mag_bits)
Generates a finite random number whose mantissa has precision at most
<tt>bits</tt> and whose exponent has at most <tt>mag_bits</tt> bits. The
values are distributed non-uniformly: special bit patterns are generated
with high probability in order to allow the test code to exercise corner
cases.
void fmpr_randtest_not_zero(fmpr_t x, flint_rand_t state, long bits,
long mag_bits)
Identical to <tt>fmpr_randtest</tt>, except that zero is never produced
as an output.
void fmpr_randtest_special(fmpr_t x, flint_rand_t state, long bits,
long mag_bits)
Indentical to <tt>fmpr_randtest</tt>, except that the output occasionally
is set to an infinity or NaN.
*******************************************************************************
Conversions
*******************************************************************************
int fmpr_get_mpfr(mpfr_t x, const fmpr_t y, mpfr_rnd_t rnd)
Sets the MPFR variable <tt>x</tt> to the value of <tt>y</tt>. If the
precision of <tt>x</tt> is too small to allow <tt>y</tt> to be represented
exactly, it is rounded in the specified MPFR rounding mode.
The return value indicates the direction of rounding,
following the standard convention of the MPFR library.
void fmpr_set_mpfr(fmpr_t x, const mpfr_t y)
Sets <tt>x</tt> to the exact value of the MPFR variable <tt>y</tt>.
void fmpr_set_ui(fmpr_t x, ulong c)
void fmpr_set_si(fmpr_t x, long c)
void fmpr_set_fmpz(fmpr_t x, const fmpz_t c)
Sets <tt>x</tt> exactly to the integer <tt>c</tt>.
void fmpr_get_fmpq(fmpq_t y, const fmpr_t x);
Sets y to the exact value of x. The result is undefined
if x is not a finite fraction.
long fmpr_set_fmpq(fmpr_t x, const fmpq_t y, long prec, fmpr_rnd_t rnd);
Sets x to the value of y, rounded according to prec and rnd.
void fmpr_set_fmpz_2exp(fmpr_t x, const fmpz_t man, const fmpz_t exp);
void fmpr_set_si_2exp_si(fmpr_t x, long man, long exp)
void fmpr_set_ui_2exp_si(fmpr_t x, ulong man, long exp)
Sets x to $\mathrm{man} \times 2^{\mathrm{exp}}$.
void fmpr_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const fmpr_t x);
Sets man and exp to the unique integers such that
$x = \mathrm{man} \times 2^{\mathrm{exp}}$ and man is odd,
provided that x is a nonzero finite fraction.
If x is zero, both man and exp are set to zero. If x is infinite or NaN,
the result is undefined.
int fmpr_get_fmpz_fixed_fmpz(fmpz_t y, const fmpr_t x, const fmpz_t e)
int fmpr_get_fmpz_fixed_si(fmpz_t y, const fmpr_t x, long e)
Converts x to a mantissa with predetermined exponent, i.e. computes
an integer y such that $y \times 2^e \approx x$, truncating if necessary.
Returns 0 if exact and 1 if truncation occurred.
*******************************************************************************
Input and output
*******************************************************************************
void fmpr_print(const fmpr_t x)
Prints the mantissa and exponent of x as integers, precisely showing
the internal representation.
void fmpr_printd(const fmpr_t x, long digits)
Prints x as a decimal floating-point number, rounding to the specified
number of digits. This function is currently implemented using MPFR,
and does not support large exponents.
*******************************************************************************
Arithmetic
*******************************************************************************
void fmpr_neg(fmpr_t y, const fmpr_t x)
Sets y to the negation of x.
long fmpr_neg_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets y to the negation of x, rounding the result.
void fmpr_abs(fmpr_t y, const fmpr_t x)
Sets y to the absolute value of x.
long fmpr_add(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_add_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_add_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_add_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = x + y$, rounded according to prec and rnd. The precision
can be FMPR_PREC_EXACT to perform an exact addition, provided that the
result fits in memory.
long _fmpr_add_eps(fmpr_t z, const fmpr_t x, int sign, long prec, fmpr_rnd_t rnd)
Sets <tt>z</tt> to the value that results by adding an infinitesimal quantity
of the given sign to <tt>x</tt>, and rounding. The result is undefined
if <tt>x</tt> is zero.
long fmpr_sub(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_sub_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_sub_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_sub_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = x - y$, rounded according to prec and rnd. The precision
can be FMPR_PREC_EXACT to perform an exact addition, provided that the
result fits in memory.
long fmpr_mul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_mul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_mul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_mul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = x \times y$, rounded according to prec and rnd. The precision
can be FMPR_PREC_EXACT to perform an exact multiplication, provided that the
result fits in memory.
void fmpr_mul_2exp_si(fmpr_t y, const fmpr_t x, long e)
void fmpr_mul_2exp_fmpz(fmpr_t y, const fmpr_t x, const fmpz_t e)
Sets y to x multiplied by $2^e$ without rounding.
long fmpr_div(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_div_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_ui_div(fmpr_t z, ulong x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_div_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_si_div(fmpr_t z, long x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_div_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
long fmpr_fmpz_div(fmpr_t z, const fmpz_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_fmpz_div_fmpz(fmpr_t z, const fmpz_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = x / y$, rounded according to prec and rnd. If $y$ is zero,
$z$ is set to NaN.
long fmpr_addmul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_addmul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_addmul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_addmul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = z + x \times y$, rounded according to prec and rnd. The
intermediate multiplication is always performed without roundoff. The
precision can be FMPR_PREC_EXACT to perform an exact addition, provided that the
result fits in memory.
long fmpr_submul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
long fmpr_submul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
long fmpr_submul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
long fmpr_submul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets $z = z - x \times y$, rounded according to prec and rnd. The
intermediate multiplication is always performed without roundoff. The
precision can be FMPR_PREC_EXACT to perform an exact subtraction, provided that the
result fits in memory.
long fmpr_sqrt(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
long fmpr_sqrt_ui(fmpr_t z, ulong x, long prec, fmpr_rnd_t rnd)
long fmpr_sqrt_fmpz(fmpr_t z, const fmpz_t x, long prec, fmpr_rnd_t rnd)
Sets $z$ to the square root of $x$, rounded according to prec and rnd.
The result is NaN if $x$ is negative.
void fmpr_pow_sloppy_fmpz(fmpr_t y, const fmpr_t b, const fmpz_t e,
long prec, fmpr_rnd_t rnd)
void fmpr_pow_sloppy_ui(fmpr_t y, const fmpr_t b, ulong e,
long prec, fmpr_rnd_t rnd)
void fmpr_pow_sloppy_si(fmpr_t y, const fmpr_t b, long e,
long prec, fmpr_rnd_t rnd)
Sets $y = b^e$, computed using without guaranteeing correct (optimal)
rounding, but guaranteeing that the result is a correct upper or lower
bound if the rounding is directional. Currently requires $b \ge 0$.
*******************************************************************************
Special functions
*******************************************************************************
long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets $z$ to $\log(x)$, rounded according to prec and rnd.
The result is NaN if $x$ is negative.
This function is currently implemented using MPFR and does not
support large exponents.
long fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets $z$ to $\log(1+x)$, rounded according to prec and rnd. This function
computes an accurate value when $x$ is small.
The result is NaN if $1+x$ is negative.
This function is currently implemented using MPFR and does not
support large exponents.
long fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets $z$ to $\exp(x)$, rounded according to prec and rnd.
This function is currently implemented using MPFR and does not
support large exponents.
long fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets $z$ to $\exp(x)-1$, rounded according to prec and rnd. This function
computes an accurate value when $x$ is small.
This function is currently implemented using MPFR and does not
support large exponents.

677
doc/fmprb.txt
View file

@ -1,677 +0,0 @@
An <tt>fmprb_t</tt> represents a ball over the real numbers,
that is, an interval $[m-r, m+r]$ where the midpoint $m$ and the
radius $r$ are (extended) real numbers and $r$ is nonnegative.
The result of an (approximate) operation done on <tt>fmprb_t</tt> variables
is a ball which contains the result of the (mathematically exact) operation
applied to any choice of points in the input balls.
In general, the output ball is not the smallest possible.
The precision parameter passed to each function roughly indicates the
precision to which calculations on the midpoint are carried out
(operations on the radius are always done using a fixed, small
precision.)
For arithmetic operations, the precision parameter currently simply
specifies the precision of the corresponding <tt>fmpr</tt> operation.
In the future, the arithmetic might be made faster by incorporating
sloppy rounding (typically equivalent to a loss of 1-2 bits of effective
working precision) when the result is known to be inexact (while still
propagating errors rigorously, of course).
Arithmetic operations done on exact input with exactly
representable output are always guaranteed to produce exact output.
For more complex operations, the precision parameter indicates a minimum
working precision (algorithms might allocate extra internal precision to
attempt to produce an output accurate to the requested number of bits,
especially when the required precision can be estimated easily, but this
is not generally required).
If the precision is increased and the inputs either are exact or are
computed with increased accuracy as well, the output should
converge proportionally, absent any bugs.
The general intended strategy for using ball arithmetic is to add a few
guard bits, and then repeat the calculation as necessary with an
exponentially increasing number of guard bits (Ziv's strategy) until the
result is exact
enough for one's purposes (typically the first attempt will be successful).
There are some caveats: in general, ball arithmetic only makes
sense for (Lipschitz) continuous functions, and
trying to approximate functions close to singularities might result in
slow convergence, or failure to converge.
<b>Warning</b>: some methods for transcendental functions and constants
currently perform the error propagation in a non-rigorous way, due to the
implementation being incomplete (in some cases, a
rigorous error bound for the algorithm or function might not
be known at all).
This should be indicated in the documentation for each function.
*******************************************************************************
Types, macros and constants
*******************************************************************************
fmprb_struct
fmprb_t
An <tt>fmprb_struct</tt> consists of a pair of <tt>fmpr_struct</tt>:s.
An <tt>fmprb_t</tt> is defined as an array of length one of type
<tt>fmprb_struct</tt>, permitting an <tt>fmprb_t</tt> to be passed by
reference.
FMPRB_RAD_PREC
The precision used for operations on the radius. This is small
enough to fit in a single word, currently 30 bits.
fmprb_midref(x)
Macro returning a pointer to the midpoint of x as an <tt>fmpr_t</tt>.
fmprb_radref(x)
Macro returning a pointer to the radius of x as an <tt>fmpr_t</tt>.
*******************************************************************************
Memory management
*******************************************************************************
void fmprb_init(fmprb_t x)
Initializes the variable x for use. Its midpoint and radius are both
set to zero.
void fmprb_clear(fmprb_t x)
Clears the variable x, freeing or recycling its allocated memory.
fmprb_struct * _fmprb_vec_init(long n)
Returns a pointer to an array of n initialized fmprb_struct:s.
void _fmprb_vec_clear(fmprb_struct * v, long n)
Clears an array of n initialized fmprb_struct:s.
*******************************************************************************
Basic manipulation
*******************************************************************************
int fmprb_is_exact(const fmprb_t x)
Returns nonzero iff the radius of x is zero.
int fmprb_equal(const fmprb_t x, const fmprb_t y)
Returns nonzero iff x and y are equal as balls, i.e. have both the
same midpoint and radius.
void fmprb_zero(fmprb_t x)
Sets x to zero.
fmprb_is_zero(const fmprb_t x)
Returns nonzero iff the midpoint and radius of x are both zero.
void fmprb_set(fmprb_t y, const fmprb_t x)
Sets y to a copy of x.
void fmprb_set_round(fmprb_t y, const fmprb_t x, long prec)
Sets y to a copy of x, rounded to prec bits.
fmprb_neg(fmprb_t y, const fmprb_t x)
Sets y to the negation of x.
fmprb_abs(fmprb_t y, const fmprb_t x)
Sets y to the absolute value of x. No attempt is made to improve the
interval represented by x if it contains zero.
void fmprb_set_fmpr(fmprb_t y, const fmpr_t x)
void fmprb_set_si(fmprb_t y, long x)
void fmprb_set_ui(fmprb_t y, ulong x)
void fmprb_set_fmpz(fmprb_t y, const fmpz_t x)
Sets y exactly to x.
void fmprb_set_fmpq(fmprb_t y, const fmpq_t x, long prec)
Sets y to the rational number x, rounded to prec bits.
int fmprb_is_one(const fmprb_t x)
Returns nonzero iff x is exactly 1.
void fmprb_one(fmprb_t x)
Sets x to the exact integer 1.
*******************************************************************************
Input and output
*******************************************************************************
void fmprb_print(const fmprb_t x)
Prints the internal representation of x.
void fmprb_printd(const fmprb_t x, long digits)
Prints x in decimal. The printed value of the radius is not adjusted
to compensate for the fact that the binary-to-decimal conversion
of both the midpoint and the radius introduces additional error.
*******************************************************************************
Random number generation
*******************************************************************************
void fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
Generates a random ball. The midpoint and radius will both be finite.
void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state,
const fmprb_t x, long bits)
Sets q to a random rational number from the interval represented by x.
A denominator is chosen by multiplying the binary denominator of x
by a random integer up to size bits.
The outcome is undefined if the midpoint or radius of x is non-finite,
or if the exponent of the midpoint or radius is so large or small
that representing the endpoints as exact rational numbers would
cause overflows.
*******************************************************************************
Precision and comparisons
*******************************************************************************
void fmprb_add_error_fmpr(fmprb_t x, const fmpr_t err)
Adds err, which is assumed to be nonnegative, to the radius of x.
void fmprb_add_error_2exp_si(fmprb_t x, long e)
Adds $2^e$ to the radius of x.
void fmprb_add_error(fmprb_t x, const fmprb_t err)
Adds the supremum of err, which is assumed to be nonnegative, to the
radius of x.
int fmprb_contains_fmpr(const fmprb_t x, const fmpr_t y)
int fmprb_contains_fmpq(const fmprb_t x, const fmpq_t y)
int fmprb_contains_fmpz(const fmprb_t x, const fmpz_t y)
int fmprb_contains_mpfr(const fmprb_t x, const mpfr_t y)
int fmprb_contains_zero(const fmprb_t x)
Returns nonzero iff the given number is contained in the interval
represented by x.
void fmprb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const fmprb_t x)
Computes the exact interval represented by x, in the form of an integer
interval multiplied by a power of two, i.e.
$x = [a, b] \times 2^{\mathrm{exp}}$.
The outcome is undefined if the midpoint or radius of x is non-finite,
or if the difference in magnitude between the midpoint and radius
is so large that representing the endpoints exactly would cause overflows.
long fmprb_rel_error_bits(const fmprb_t x)
Returns the effective relative error of self measured in bits, defined as
the difference between the position of the top bit in the radius
and the top bit in the midpoint, plus one.
The result is clamped between plus/minus FMPR_PREC_EXACT.
long fmprb_rel_accuracy_bits(const fmprb_t x)
Returns the effective relative accuracy of x measured in bits,
equal to the negative of the return value from fmprb_rel_error_bits.
*******************************************************************************
Arithmetic
*******************************************************************************
void fmprb_add(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_add_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_add_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_add_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
void fmprb_add_fmpr(fmprb_t z, const fmprb_t x, const fmpr_t y, long prec)
Sets $z = x + y$, rounded to prec bits. The precision can be
FMPR_PREC_EXACT provided that the result fits in memory.
void fmprb_sub(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_sub_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_sub_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_sub_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets $z = x - y$, rounded to prec bits. The precision can be
FMPR_PREC_EXACT provided that the result fits in memory.
void fmprb_mul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_mul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_mul_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_mul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets $z = x \times y$, rounded to prec bits. The precision can be
FMPR_PREC_EXACT provided that the result fits in memory.
void fmprb_mul_2exp_si(fmprb_t y, const fmprb_t x, long e)
Sets $y$ to $x$ multiplied by $2^e$.
void fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_div_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_div_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_div_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
void fmprb_fmpz_div_fmpz(fmprb_t y, const fmpz_t num, const fmpz_t den, long prec)
void fmprb_ui_div(fmprb_t z, ulong x, const fmprb_t y, long prec);
Sets $z = x / y$, rounded to prec bits. If $y$ contains zero, $z$ is
set to $0 \pm \infty$. Otherwise, error propagation uses the rule
$$\left| \frac{x}{y} - \frac{x+\xi_1 a}{y+\xi_2 b} \right| =
\left|\frac{x \xi_2 b - y \xi_1 a}{y (y+\xi_2 b)}\right| \le
\frac{|xb|+|ya|}{|y| (|y|-b)}$$
where $-1 \le \xi_1, \xi_2 \le 1$, and
where the triangle inequality has been applied to the numerator and
the reverse triangle inequality has been applied to the denominator.
void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec);
Sets $y = x / (2^n - 1)$, rounded to prec bits.
void fmprb_addmul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_addmul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_addmul_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_addmul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets $z = z + x \times y$, rounded to prec bits. The precision can be
FMPR_PREC_EXACT provided that the result fits in memory.
void fmprb_submul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
void fmprb_submul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
void fmprb_submul_si(fmprb_t z, const fmprb_t x, long y, long prec)
void fmprb_submul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets $z = z - x \times y$, rounded to prec bits. The precision can be
FMPR_PREC_EXACT provided that the result fits in memory.
void fmprb_sqrt(fmprb_t z, const fmprb_t x, long prec)
void fmprb_sqrt_ui(fmprb_t z, ulong x, long prec)
void fmprb_sqrt_fmpz(fmprb_t z, const fmpz_t x, long prec)
Sets $z$ to the square root of $x$, rounded to prec bits.
Error propagation is done using the following rule:
assuming $m > r \ge 0$, the error is largest at $m - r$, and we have
$\sqrt{m} - \sqrt{m-r} \le r / (2 \sqrt{m-r})$.
void fmprb_pow_fmpz(fmprb_t y, const fmprb_t b, const fmpz_t e, long prec)
void fmprb_pow_ui(fmprb_t y, const fmprb_t b, ulong e, long prec)
void fmprb_ui_pow_ui(fmprb_t y, ulong b, ulong e, long prec)
void fmprb_si_pow_ui(fmprb_t y, long b, ulong e, long prec)
Sets $y = b^e$ using binary exponentiation. Provided that $b$ and $e$
are small enough and the exponent is positive, the exact power can be
computed using FMPR_PREC_EXACT.
*******************************************************************************
Special functions
*******************************************************************************
void fmprb_log(fmprb_t z, const fmprb_t x, long prec)
void fmprb_log_ui(fmprb_t z, ulong x, long prec)
void fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec)
Sets $z = \log(x)$. Error propagation is done using the following rule:
assuming $m > r \ge 0$, the error is largest at $m - r$, and we have
$\log(m) - \log(m-r) = \log(1 + r/(m-r))$. The last expression is
calculated accurately for small radii via <tt>fmpr_log1p</tt>.
An input containing zero currently raises an exception.
void fmprb_exp(fmprb_t z, const fmprb_t x, long prec)
Sets $z = \exp(x)$. Error propagation is done using the following rule:
the error is largest at $m + r$, and we have
$\exp(m+r) - \exp(m) = \exp(m) (\exp(r)-1)$.
The last expression is calculated accurately for small radii
via <tt>fmpr_expm1</tt>.
void fmprb_sin(fmprb_t s, const fmprb_t x, long prec)
void fmprb_cos(fmprb_t c, const fmprb_t x, long prec)
void fmprb_sin_cos(fmprb_t s, fmprb_t c, const fmprb_t x, long prec)
Sets $s = \sin x$, $c = \cos x$. Error propagation uses the rule
$|\sin(m \pm r) - \sin(m)| \le r$ (this could be tightened to
$\min(r,2)$).
void fmprb_atan(fmprb_t z, const fmprb_t x, long prec)
Sets $z = \tan^{-1} x$. Letting $d = \max(0, |m| - r)$,
the propagated error is bounded by $r / (1 + d^2)$
(this could be tightened).
void fmprb_fac_ui(fmprb_t x, ulong n, long prec)
Sets x to n factorial, computed using binary splitting. Provided that
$n$ is small enough, the exact factorial can be computed using
FMPR_PREC_EXACT.
void fmprb_rfac_ui_bsplit(fmprb_t y, const fmprb_t x, ulong n, long prec)
Sets x to the rising factorial $x (x+1) (x+2) \cdots (x+n-1)$,
computed using binary splitting.
The basecase processes eight factors at a time using the formula
$x(x+1)\cdots(x+7) = (28 + 98x + 63x^2 + 14x^3 + x^4)^2 - 16 (7+2x)^2$,
replacing 7 full-precision multiplications with 4 squarings and
1 multiplication [1].
Empirically, this is about twice as fast at high precision.
Numerical stability is slightly worse.
[1] R. Crandall and C. Pomerance, "Prime Numbers: A Computational
Perspective", second edition, Springer (2005), p. 316.
void fmprb_rfac_ui_multipoint(fmprb_t y, const fmprb_t x, ulong n, long prec)
Sets x to the rising factorial $x (x+1) (x+2) \cdots (x+n-1)$,
computed using fast multipoint evaluation. This only requires
$O(n^{1/2+\varepsilon})$ multiplications, but has high overhead
and poor numerical stability (adding $O(n)$ guard bits to the input
might be necessary to achieve full accuracy). It can be expected to
be faster than the binary splitting algorithm if the input is a
full-precision number, the precision is at least 100000 bits,
and $n$ is of the same order of magnitude as (perhaps slightly
smaller than) the number of bits.
void fmprb_bin_ui(fmprb_t x, const fmprb_t n, ulong k, long prec)
void fmprb_bin_uiui(fmprb_t x, ulong n, ulong k, long prec)
Sets x to the binomial coefficient ${n \choose k}$, computed using
binary splitting. Provided that $n$ and $k$ are small enough, an exact
binomial coefficient can be computed using FMPR_PREC_EXACT.
void fmprb_fib_fmpz(fmprb_t f, const fmpz_t n, long prec)
void fmprb_fib_ui(fmprb_t f, ulong n, long prec)
Sets x to the Fibonacci number $F_n$. Uses the binary squaring
algorithm described in D. Takahashi, "A fast algorithm for
computing large Fibonacci numbers",
Information Processing Letters 75 (2000) 243-246
Provided that $n$ is small enough, an exact Fibonacci number can be
computed using FMPR_PREC_EXACT.
void fmprb_const_pi_chudnovsky(fmprb_t x, long prec)
Sets x to $\pi$, computed using the Chudnovsky algorithm.
Letting $A = 13591409$, $B = 545140134$, $C = 640320$,
we have $\pi \approx 1 / s_N$ where
$$s_N = 12 \sum_{k=0}^N \frac{(-1)^k (6k)! (A+Bk)}
{(3k)! (k!)^3 C^{3k+3/2}}$$
The implementation computes an approximation for the
algebraic number $1/s_N$ using binary splitting, bounding
the rounding error automatically.
The hypergeometric term ratio is asymptotically
$R = C^3 / (2^6 \times 3^3) \approx 1.5 \times 10^{14}$, and in fact we have
$|\pi - 1/s_N| < 1/R^N$ (with a more detailed calculation, the truncation
error could be bounded closer to $1/R^{N+1}$).
void fmprb_const_pi(fmprb_t x, long prec)
Sets x to $\pi$. The value is cached for repeated use.
void fmprb_const_log_sqrt2pi(fmprb_t x, long prec)
Sets x to $\log \sqrt{2 \pi}$. The value is cached for repeated use.
void fmprb_const_euler_brent_mcmillan(fmprb_t res, long prec)
Sets x to Euler's constant $\gamma$, computed using the second
Bessel function formula of Brent and McMillan.
Brent and McMillan conjectured that the error depending
on the internal parameter $n$ is of order $O(e^{-8n})$. Brent has
recently proved that this bound is correct, but without determining
an explicit big-O factor.
[1] R. P. Brent and E. M. McMillan,
"Some new algorithms for high-precision computation of Euler's constant",
Mathematics of Computation 34 (1980), 305-312.</li>
[2] R. P. Brent, "Ramanujan and Euler's Constant",
http://wwwmaths.anu.edu.au/~brent/pd/Euler_CARMA_10.pdf</li>
[3] The MPFR team (2012), "MPFR Algorithms",
http://www.mpfr.org/algo.html</li>
void fmprb_const_zeta3_bsplit(fmprb_t x, long prec)
Sets x to Apery's constant $\zeta(3)$, computed by applying binary
splitting to a hypergeometric series.
void fmprb_zeta_ui_asymp(fmprb_t z, ulong s, long prec)
Assuming $s \ge 2$, approximates $\zeta(s)$ by $1 + 2^{-s}$ along with
a correct error bound. We use the following bounds: for $s > b$,
$\zeta(s) - 1 < 2^{-b}$, and generally,
$\zeta(s) - (1 + 2^{-s}) < 2^{2-\lfloor 3 s/2 \rfloor}$.
void fmprb_zeta_ui_euler_product(fmprb_t z, ulong s, long prec)
Computes $\zeta(s)$ using the Euler product. This is fast only if s
is large compared to the precision.
Writing $P(a,b) = \prod_{a \le p \le b} (1 - p^{-s})$, we have
$1/\zeta(s) = P(a,M) P(M+1,\infty)$.
To bound the error caused by truncating
the product at $M$, we write $P(M+1,\infty) = 1 - \epsilon(s,M)$.
Since $0 < P(a,M) \le 1$, the absolute error for $\zeta(s)$ is
bounded by $\epsilon(s,M)$.
According to the analysis in [1], it holds for all $s \ge 6$ and $M \ge 1$
that $1/P(M+1,\infty) - 1 \le f(s,M) \equiv 2 M^{1-s} / (s/2 - 1)$.
Thus, we have $1/(1-\epsilon(s,M)) - 1 \le f(s,M)$, and expanding
the geometric series allows us to conclude that
$\epsilon(M) \le f(s,M)$.
[1] S. Fillebrown, "Faster Computation of Bernoulli Numbers",
Journal of Algorithms 13, 431-445 (1992)
void fmprb_zeta_ui_bernoulli(fmprb_t x, ulong n, long prec)
Computes $\zeta(n)$ for even $n$ via the corresponding Bernoulli number,
which is generated using FLINT.
void fmprb_zeta_ui_vec_borwein(fmprb_struct * z, ulong start, long num,
ulong step, long prec)
Evaluates $\zeta(s)$ at $\mathrm{num}$ consecutive integers $s$ beginning
with $\mathrm{start}$ and proceeding in increments of $\mathrm{step}$.
Uses Borwein's formula [1], implemented to support fast multi-evaluation
(but also works well for a single $s$).
Requires $\mathrm{start} \ge 2$. For efficiency, the largest $s$
should be at most about as
large as $\mathrm{prec}$. Arguments approaching LONG_MAX will cause
overflows.
One should therefore only use this function for s up to about prec, and
then switch to the Euler product.
The algorithm for single $s$ is basically identical to the one used in MPFR
(see [2] for a detailed description).
In particular, we evaluate the sum backwards to avoid storing more than one
$d_k$ coefficient, and use integer arithmetic throughout since it
is convenient and the terms turn out to be slightly larger than
$2^\mathrm{prec}$.
The only numerical error in the main loop comes from the division by $k^s$,
which adds less than 1 unit of error per term.
For fast multi-evaluation, we repeatedly divide by $k^{\mathrm{step}}$.
Each division reduces the input error and adds at most 1 unit of
additional rounding error, so by induction, the error per term
is always smaller than 2 units.
References:
[1] P. Borwein, "An Efficient Algorithm for the Riemann Zeta Function",
Constructive experimental and nonlinear analysis,
CMS Conference Proc. 27 (2000), 29-34
http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
[2] The MPFR team (2012), "MPFR Algorithms", http://www.mpfr.org/algo.html
[3] X. Gourdon and P. Sebah (2003),
"Numerical evaluation of the Riemann Zeta-function"
http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf
void fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec)
Computes $\zeta(s)$ for arbitrary $s \ge 2$ using a binary splitting
implementation of Borwein's formula. The algorithm has quasilinear
complexity with respect to the precision.
void fmprb_zeta_ui(fmprb_t x, ulong s, long prec)
Computes $\zeta(s)$ for nonnegative integer $s \ne 1$, automatically
choosing an appropriate algorithm.
void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec)
void fmprb_zeta_ui_vec_even(fmprb_struct * x, ulong start, long num, long prec)
void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec)
Computes $\zeta(s)$ at num consecutive integers (respectively num
even or num odd integers) beginning with $s = \mathrm{start} \ge 2$,
automatically choosing an appropriate algorithm.
void fmprb_gamma_fmpq_karatsuba(fmprb_struct * v, const fmpq_t a,
long num, long prec)
Uses Karatsuba's algorithm [1] to compute num coefficients in the
Taylor series of $\Gamma(a+x)$ for rational $0 < a &le; 1$, i.e.
computes $\Gamma(a), \Gamma'(a) ... \Gamma^{(\mathrm{num}-1)}(a) / (\mathrm{num}-1)!$
This algorithm is most efficient at high precision, for num much smaller
than the number of bits, and with small denominators of $a$.
In particular, with num = 1, this algorithm computes $\Gamma(a)$
efficiently for small rational $a$.
Let $s = \max(2, \mathrm{num}-1)$. With parameters $r$ and $n$ chosen
such that $r \ge n$ and $n \ge 2 s \log 2 s$, Karatsuba shows that
$$\Gamma^{(j)}(a) = \sum_{k=0}^r \frac{(-1)^k}{k!}
\frac{n^{k+a}}{k+a}
\sum_{m=0}^j (-1)^m \frac{j! \, \log^{j-m} n}{(j-m)! (k+a)^m}
+ \theta_j$$
where $$|\theta_j| \le \frac{5}{3} e^{-n} \log^s n +
\left(\frac{e}{r+2}\right)^{r+2} (1 + n^{r+2} \log^s n).$$
We choose the parameters $n$ and $r$ heuristically to be nearly optimal,
and then evaluate the above formula to bound $\theta_j$ rigorously.
Karatsuba claims that choosing $r \ge 3n$ gives
$|\theta_j| \le 2^{-n-1}$. This is, unfortunately, incorrect.
Setting $r = n \alpha$ and expanding the error term around $n = \infty$,
one finds that $\alpha$ asymptotically should be
$1/W(1/e) \approx 3.59112147666862$ where $W(x)$ is the Lambert W-function.
We also optimize the selection of $n$ by choosing
$n \approx b \log 2$ where $b$ is the desired number of bits, rather
than $n \approx b$, and round $n$ so that it has a short binary expansion
(this gives smaller numbers in the binary splitting stage).
Finally, if $s$ is small, we perform binary splitting to a working
precision of about $2.2$ times the target precision rather than exactly.
This factor was tested to give full accuracy up to at least one million
digits when $s \approx 1$. A more careful analysis should
be done here so that a working precision is selected which always is
sufficient and also nearly optimal.
[1] E. A. Karatsuba, "Fast evaluation of the Hurwitz zeta function and
Dirichlet L-series", Problems of Information Transmission,
Vol. 34, No. 4, 1998.
void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec)
Sets $y = \log \Gamma(x)$, assuming that $x > 0$.
For large $x$, uses Stirling's expansion
$$\log \Gamma(x) = \left(x-\frac{1}{2}\right)\log x - x +
\frac{\ln {2 \pi}}{2} + \sum_{k=1}^{n-1} t_k + R(n,x)
$$
where $$t_k = \frac{B_{2k}}{2k(2k-1)x^{2k-1}}$$ and $|R(n,x)| < t_n$.
If $x$ is too small for the asymptotic expansion to give sufficient
accuracy directly, we translate to $x + r$
using the formula $\log \Gamma(x) = \log \Gamma(x+r) -
\log(x (x+1) (x+2) \cdots (x+r-1))$.
To obtain a remainder smaller than $2^{-b}$, we must choose an $r$ such
that $x + r > \beta b$, where $\beta > \log(2) / (2 \pi) \approx 0.11$.
We use a slightly larger factor $\beta \approx 0.2$ to more closely
balance $n$ and $r$. A much larger $\beta$ (e.g. $\beta = 1$) could be
used to reduce the number of Bernoulli numbers that have to be
precomputed, at the expense of slower repeated evaluation.

203
doc/fmprb_mat.txt
View file

@ -1,203 +0,0 @@
An <tt>fmprb_mat_t</tt> represents a dense matrix over the real numbers,
implemented as an array of entries of type <tt>fmprb_struct</tt>.
The dimension (number of rows and columns) of a matrix is fixed at
initialization, and the user must ensure that inputs and outputs to
an operation have compatible dimensions. The number of rows or columns
in a matrix can be zero.
*******************************************************************************
Types, macros and constants
*******************************************************************************
fmprb_mat_struct
fmprb_mat_t
Contains a pointer to a flat array of the entries (entries), an array of
pointers to the start of each row (rows), and the number of rows (r)
and columns (c).
An <tt>fmprb_mat_t</tt> is defined as an array of length one of type
<tt>fmprb_mat_struct</tt>, permitting an <tt>fmprb_mat_t</tt> to
be passed by reference.
fmprb_mat_entry(mat, i, j)
Macro giving a pointer to the entry at row i and column j.
fmprb_mat_nrows(mat)
Returns the number of rows of the matrix.
fmprb_mat_ncols(mat)
Returns the number of columns of the matrix.
*******************************************************************************
Memory management
*******************************************************************************
void fmprb_mat_init(fmprb_mat_t mat, long r, long c)
Initializes the matrix, setting it to the zero matrix with the given
number of rows and columns.
void fmprb_mat_clear(fmprb_mat_t mat)
Clears the matrix, deallocating all entries.
*******************************************************************************
Conversions
*******************************************************************************
void fmprb_mat_set(fmprb_mat_t dest, const fmprb_mat_t src)
void fmprb_mat_set_fmpz_mat(fmprb_mat_t dest, const fmpz_mat_t src)
void fmprb_mat_set_fmpq_mat(fmprb_mat_t dest, const fmpq_mat_t src, long prec)
Sets dest to src. The operands must have identical dimensions.
*******************************************************************************
Input and output
*******************************************************************************
void fmprb_mat_printd(const fmprb_mat_t mat, long digits)
Prints each entry in the matrix with the specified number of decimal digits.
*******************************************************************************
Comparisons
*******************************************************************************
int fmprb_mat_equal(const fmprb_mat_t mat1, const fmprb_mat_t mat2)
Returns nonzero iff the matrices have the same dimensions
and identical entries.
int fmprb_mat_contains_fmpq_mat(const fmprb_mat_t mat1, const fmpq_mat_t mat2)
Returns nonzero iff the matrices have the same dimensions and each entry
in mat2 is contained in the corresponding entry in mat1.
*******************************************************************************
Special matrices
*******************************************************************************
void fmprb_mat_zero(fmprb_mat_t mat)
Sets all entries in mat to zero.
void fmprb_mat_one(fmprb_mat_t mat)
Sets the entries on the main diagonal to ones,
and all other entries to zero.
*******************************************************************************
Arithmetic
*******************************************************************************
void fmprb_mat_neg(fmprb_mat_t dest, const fmprb_mat_t src)
Sets dest to the exact negation of src. The operands must have
the same dimensions.
void fmprb_mat_add(fmprb_mat_t res, const fmprb_mat_t mat1,
const fmprb_mat_t mat2, long prec)
Sets res to the sum of mat1 and mat2. The operands must have
the same dimensions.
void fmprb_mat_sub(fmprb_mat_t res, const fmprb_mat_t mat1,
const fmprb_mat_t mat2, long prec)
Sets res to the difference of mat1 and mat2. The operands must have
the same dimensions.
void fmprb_mat_mul(fmprb_mat_t res, const fmprb_mat_t mat1,
const fmprb_mat_t mat2, long prec)
Sets res to the matrix product of mat1 and mat2. The operands must have
compatible dimensions for matrix multiplication.
*******************************************************************************
Solving
*******************************************************************************
int fmprb_mat_lu(long * perm, fmprb_mat_t LU, const fmprb_mat_t A, long prec)
Given an $n \times n$ matrix $A$, computes an LU decomposition $PLU = A$
using Gaussian elimination with partial pivoting.
The input and output matrices can be the same, performing the
decomposition in-place.
Entry $i$ in the permutation vector perm is set to the row index in
the input matrix corresponding to row $i$ in the output matrix.
The algorithm succeeds and returns nonzero if it can find $n$ invertible
(i.e. not containing zero) pivot entries. This guarantees that the matrix
is invertible.
The algorithm fails and returns zero, leaving the entries in $P$ and $LU$
undefined, if it cannot find $n$ invertible pivot elements.
In this case, either the matrix is singular, the input matrix was
computed to insufficient precision, or the LU decomposition was
attempted at insufficient precision.
void fmprb_mat_solve_lu_precomp(fmprb_mat_t X, const long * perm,
const fmprb_mat_t LU, const fmprb_mat_t B, long prec)
Solves $AX = B$ given the precomputed nonsingular LU decomposition $A = PLU$.
The matrices $X$ and $B$ are allowed to be aliased with each other,
but $X$ is not allowed to be aliased with $LU$.
int fmprb_mat_solve(fmprb_mat_t X, const fmprb_mat_t A,
const fmprb_mat_t B, long prec)
Solves $AX = B$ where $A$ is a nonsingular $n \times n$ matrix
and $X$ and $B$ are $n \times m$ matrices, using LU decomposition.
If $m > 0$ and $A$ cannot be inverted numerically (indicating either that
$A$ is singular or that the precision is insufficient), the values in the
output matrix are left undefined and zero is returned. A nonzero return
value guarantees that $A$ is invertible and that the exact solution
matrix is contained in the output.
int fmprb_mat_inv(fmprb_mat_t X, const fmprb_mat_t A, long prec)
Sets $X = A^{-1}$ where $A$ is a square matrix, computed by solving
the system $AX = I$.
If $A$ cannot be inverted numerically (indicating either that
$A$ is singular or that the precision is insufficient), the values in the
output matrix are left undefined and zero is returned.
A nonzero return value guarantees that the matrix is invertible
and that the exact inverse is contained in the output.
void fmprb_mat_det(fmprb_t det, const fmprb_mat_t A, long prec)
Computes the determinant of the matrix, using Gaussian elimination
with partial pivoting. If at some point an invertible pivot element
cannot be found, the elimination is stopped and the magnitude of the
determinant of the remaining submatrix is bounded using
Hadamard's inequality.

431
doc/fmprb_poly.txt
View file

@ -1,431 +0,0 @@
An <tt>fmprb_poly_t</tt> represents a polynomial over the real numbers,
implemented as an array of coefficients of type <tt>fmprb_struct</tt>.
Most functions are provided in two versions: an underscore method which
operates directly on pre-allocated arrays of coefficients and generally
has some restrictions (such as requiring the lengths to be nonzero
and not supporting aliasing of the input and output arrays),
and a non-underscore method which performs automatic memory
management and handles degenerate cases.
*******************************************************************************
Types, macros and constants
*******************************************************************************
fmprb_poly_struct
fmprb_poly_t
Contains a pointer to an array of coefficients (coeffs), the used
length (length), and the allocated size of the array (alloc).
An <tt>fmprb_poly_t</tt> is defined as an array of length one of type
<tt>fmprb_poly_struct</tt>, permitting an <tt>fmprb_poly_t</tt> to
be passed by reference.
*******************************************************************************
Memory management
*******************************************************************************
void fmprb_poly_init(fmprb_poly_t poly)
Initializes the polynomial for use, setting it to the zero polynomial.
void fmprb_poly_clear(fmprb_poly_t poly)
Clears the polynomial, deallocating all coefficients and the
coefficient array.
void fmprb_poly_fit_length(fmprb_poly_t poly, long len)
Makes sures that the coefficient array of the polynomial contains at
least len initialized coefficients.
void _fmprb_poly_set_length(fmprb_poly_t poly, long len)
Directly changes the length of the polynomial, without allocating or
deallocating coefficients. The value shold not exceed the allocation length.
void _fmprb_poly_normalise(fmprb_poly_t poly)
Strips any trailing coefficients which are identical to zero.
The output
void fmprb_poly_zero(fmprb_poly_t poly)
void fmprb_poly_one(fmprb_poly_t poly)
Sets poly to the constant 0 respectively 1.
*******************************************************************************
Conversions
*******************************************************************************
void fmprb_poly_set_fmpz_poly(fmprb_poly_t poly, const fmpz_poly_t src, long prec)
void fmprb_poly_set_fmpq_poly(fmprb_poly_t poly, const fmpq_poly_t src, long prec)
Sets poly to src.
*******************************************************************************
Input and output
*******************************************************************************
void fmprb_poly_printd(const fmprb_poly_t poly, long digits)
Prints the polynomial as an array of coefficients, printing each
coefficient using fmprb_printd.
*******************************************************************************
Comparisons
*******************************************************************************
int fmprb_poly_contains_fmpq_poly(const fmprb_poly_t poly1,
const fmpq_poly_t poly2)
Returns nonzero iff poly1 contains poly2.
int fmprb_poly_equal(const fmprb_t A, const fmprb_t B)
Returns nonzero iff A and B are equal as polynomial balls, i.e. all
coefficients have equal midpoint and radius.
*******************************************************************************
Arithmetic
*******************************************************************************
void
_fmprb_poly_add(fmprb_struct * C, const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long prec)
Sets {C, max(lenA, lenB)} to the sum of {A, lenA} and {B, lenB}.
Allows aliasing of the input and output operands.
void fmprb_poly_add(fmprb_poly_t C, const fmprb_poly_t A,
const fmprb_poly_t B, long prec)
Sets C to the sum of A and B.
void _fmprb_poly_mullow(fmprb_struct * C,
const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long n, long prec)
Sets {C, n} to the product of {A, lenA} and {B, lenB}, truncated to
length n. The output is not allowed to be aliased with either of the
inputs. We require lenA &ge; lenB > 0, n > 0, lenA + lenB - 1 &ge; n.
As currently implemented, this function puts each input polynomial on
a common exponent, truncates to prec bits, and multiplies exactly over
the integers. The output error is computed by cross-multiplying the
max norms.
void fmprb_poly_mullow(fmprb_poly_t C, const fmprb_poly_t A,
const fmprb_poly_t B, long n, long prec)
Sets C to the product of A and B, truncated to length n.
void _fmprb_poly_mul(fmprb_struct * C,
const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long prec)
Sets {C, n} to the product of {A, lenA} and {B, lenB}, truncated to
length n. The output is not allowed to be aliased with either of the
inputs. We require lenA &ge; lenB > 0, n > 0.
This function currently calls _fmprb_poly_mullow.
void fmprb_poly_mul(fmprb_poly_t C, const fmprb_poly_t A,
const fmprb_poly_t B, long prec)
Sets C to the product of A and B.
void _fmprb_poly_inv_series(fmprb_struct * Qinv, const fmprb_struct * Q,
long len, long prec)
Sets {Qinv, len} to the power series inverse of {Q, len}.
Uses Newton iteration.
void fmprb_poly_inv_series(fmprb_poly_t Qinv, const fmprb_poly_t Q,
long n, long prec)
Sets Qinv to the power series inverse of Q.
void _fmprb_poly_div(fmprb_struct * Q,
const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long prec)
void _fmprb_poly_rem(fmprb_struct * R,
const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long prec)
void _fmprb_poly_divrem(fmprb_struct * Q, fmprb_struct * R,
const fmprb_struct * A, long lenA,
const fmprb_struct * B, long lenB, long prec)
void fmprb_poly_divrem(fmprb_poly_t Q, fmprb_poly_t R,
const fmprb_poly_t A, const fmprb_poly_t B, long prec)
Performs polynomial division with remainder, computing a quotient $Q$ and
a remainder $R$ such that $A = BQ + R$. The leading coefficient of $B$ must
not contain zero. The implementation reverses the inputs and performs
power series division.
void _fmprb_poly_div_root(fmprb_struct * Q, fmprb_t R, const fmprb_struct * A,
long len, const fmprb_t c, long prec)
Divides $A$ by the polynomial $x - c$, computing the quotient $Q$ as well
as the remainder $R = f(c)$.
*******************************************************************************
Product trees
*******************************************************************************
void _fmprb_poly_product_roots(fmprb_struct * poly, const fmprb_struct * xs, long n, long prec)
void fmprb_poly_product_roots(fmprb_poly_t poly, fmprb_struct * xs, long n, long prec)
Generates the polynomial $(x-x_0)(x-x_1)\cdots(x-x_{n-1})$.
fmprb_struct ** _fmprb_poly_tree_alloc(long len)
Returns an initialized data structured capable of representing a
remainder tree (product tree) of the given number of roots.
void _fmprb_poly_tree_free(fmprb_struct ** tree, long len)
Deallocates a tree structure as allocated using _fmprb_poly_tree_alloc.
void _fmprb_poly_tree_build(fmprb_struct ** tree,
const fmprb_struct * roots, long len, long prec)
Constructs a product tree from a given array of roots. The tree
structure must be pre-allocated to the specified length using
_fmprb_poly_tree_alloc.
*******************************************************************************
Composition
*******************************************************************************
void _fmprb_poly_compose_horner(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long prec)
void fmprb_poly_compose_horner(fmprb_poly_t res,
const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
void _fmprb_poly_compose_divconquer(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long prec);
void fmprb_poly_compose_divconquer(fmprb_poly_t res,
const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
void _fmprb_poly_compose(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long prec)
void fmprb_poly_compose(fmprb_poly_t res,
const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
Sets res to the composition $h(x) = f(g(x))$ where $f$ is given by
poly1 and $g$ is given by poly2, respectively using Horner's rule,
divide-and-conquer, and an automatic choice between the two algorithms.
The underscore methods do not support aliasing of the output
with either input polynomial.
void _fmprb_poly_compose_series_horner(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long n, long prec)
void fmprb_poly_compose_series_horner(fmprb_poly_t res,
const fmprb_poly_t poly1,
const fmprb_poly_t poly2, long n, long prec)
void _fmprb_poly_compose_series_brent_kung(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long n, long prec)
void fmprb_poly_compose_series_brent_kung(fmprb_poly_t res,
const fmprb_poly_t poly1,
const fmprb_poly_t poly2, long n, long prec)
void _fmprb_poly_compose_series(fmprb_struct * res,
const fmprb_struct * poly1, long len1,
const fmprb_struct * poly2, long len2, long n, long prec)
void fmprb_poly_compose_series(fmprb_poly_t res,
const fmprb_poly_t poly1,
const fmprb_poly_t poly2, long n, long prec)
Sets res to the power series composition $h(x) = f(g(x))$ truncated
to order $O(x^n)$ where $f$ is given by poly1 and $g$ is given by poly2,
respectively using Horner's rule, the Brent-Kung baby step-giant step
algorithm, and an automatic choice between the two algorithms.
We require that the constant term in $g(x)$ is exactly zero.
The underscore methods do not support aliasing of the output
with either input polynomial.
*******************************************************************************
Evaluation and interpolation
*******************************************************************************
void _fmprb_poly_evaluate(fmprb_t res, const fmprb_struct * f, long len,
const fmprb_t a, long prec)
void fmprb_poly_evaluate(fmprb_t res, const fmprb_poly_t f, const fmprb_t a, long prec)
Sets res to $f(a)$, evaluated using Horner's rule.
void _fmprb_poly_evaluate_vec_iter(fmprb_struct * ys, const fmprb_struct * poly, long plen,
const fmprb_struct * xs, long n, long prec)
void fmprb_poly_evaluate_vec_iter(fmprb_struct * ys,
const fmprb_poly_t poly, const fmprb_struct * xs, long n, long prec)
Evaluates the polynomial simultaneously at n given points, calling
_fmprb_poly_evaluate repeatedly.
void _fmprb_poly_evaluate_vec_fast_precomp(fmprb_struct * vs, const fmprb_struct * poly,
long plen, fmprb_struct ** tree, long len, long prec)
void _fmprb_poly_evaluate_vec_fast(fmprb_struct * ys, const fmprb_struct * poly, long plen,
const fmprb_struct * xs, long n, long prec)
void fmprb_poly_evaluate_vec_fast(fmprb_struct * ys,
const fmprb_poly_t poly, const fmprb_struct * xs, long n, long prec)
Evaluates the polynomial simultaneously at n given points, using
fast multipoint evaluation.
void _fmprb_poly_interpolate_newton(fmprb_struct * poly, const fmprb_struct * xs,
const fmprb_struct * ys, long n, long prec)
void fmprb_poly_interpolate_newton(fmprb_poly_t poly,
const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most n that interpolates
the given x and y values. This implementation first interpolates in the
Newton basis and then converts back to the monomial basis.
void _fmprb_poly_interpolate_barycentric(fmprb_struct * poly,
const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
void fmprb_poly_interpolate_barycentric(fmprb_poly_t poly,
const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most n that interpolates
the given x and y values. This implementation uses the barycentric
form of Lagrange interpolation.
void _fmprb_poly_interpolation_weights(fmprb_struct * w,
fmprb_struct ** tree, long len, long prec)
void _fmprb_poly_interpolate_fast_precomp(fmprb_struct * poly,
const fmprb_struct * ys, fmprb_struct ** tree, const fmprb_struct * weights,
long len, long prec)
void _fmprb_poly_interpolate_fast(fmprb_struct * poly,
const fmprb_struct * xs, const fmprb_struct * ys, long len, long prec)
void fmprb_poly_interpolate_fast(fmprb_poly_t poly,
const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most n that interpolates
the given x and y values, using fast Lagrange interpolation.
The precomp function takes a precomputed product tree over the
x values and a vector of interpolation weights as additional inputs.
*******************************************************************************
Differentiation
*******************************************************************************
void _fmprb_poly_derivative(fmprb_struct * res, const fmprb_struct * poly,
long len, long prec)
Sets {res, len - 1} to the derivative of {poly, len}.
Allows aliasing of the input and output.
void fmprb_poly_derivative(fmprb_poly_t res, const fmprb_poly_t poly, long prec)
Sets res to the derivative of poly.
void _fmprb_poly_integral(fmprb_struct * res, const fmprb_struct * poly,
long len, long prec)
Sets {res, len} to the integral of {poly, len - 1}.
Allows aliasing of the input and output.
void fmprb_poly_integral(fmprb_poly_t res, const fmprb_poly_t poly, long prec)
Sets res to the integral of poly.
*******************************************************************************
Special functions
*******************************************************************************
void _fmprb_poly_log_series(fmprb_struct * f, fmprb_struct * h,
long n, long prec)
void fmprb_poly_log_series(fmprb_poly_t f, const fmprb_poly_t h,
long n, long prec)
Sets $f$ to the power series logarithm of $h$, truncated to length $n$.
Uses the formula $\log f = \int f' / f$, adding the logarithm of the
constant term in $h$ as the constant of integration.
The underscore method does not support aliasing of the input and output
arrays.
void _fmprb_poly_exp_series_basecase(fmprb_struct * f,
const fmprb_struct * h, long hlen, long n, long prec)
void fmprb_poly_exp_series_basecase(fmprb_poly_t f, const fmprb_poly_t h,
long n, long prec)
void _fmprb_poly_exp_series(fmprb_struct * f, const fmprb_struct * h, long hlen, long n, long prec)
void fmprb_poly_exp_series(fmprb_poly_t f, const fmprb_poly_t h, long n, long prec)
Sets $f$ to the power series exponential of $h$, truncated to length $n$.
The basecase version uses a simple recurrence for the coefficients,
requiring $O(nm)$ operations where $m$ is the length of $h$.
The main implementation uses Newton iteration, starting from a small
number of terms given by the basecase algorithm. The complexity
is $O(M(n))$. Redundant operations in the Newton iteration are
avoided by using the scheme described in [1].
The underscore methods support aliasing and allow the input to be
shorter than the output, but require the lengths to be nonzero.
[1] G. Hanrot and P. Zimmermann, "Newton Iteration Revisited",
http://www.loria.fr/~zimmerma/papers/fastnewton.ps.gz (2004 preprint)
void fmprb_poly_log_gamma_series(fmprb_poly_t f, long n, long prec)
Sets $f$ to the series expansion of $\log(\Gamma(1-x))$, truncated to
length $n$.

24
doc/history.txt
View file

@ -1,24 +0,0 @@
<ul>
<li>2012-09-29 - version 0.2
<ul>
<li>code for computing the gamma function (Karatsuba, Stirling's series)</li>
<li>rising factorials</li>
<li>fast exp_series using Newton iteration</li>
<li>improved multiplication of small polynomials by using classical multiplication</li>
<li>implemented error propagation for square roots</li>
<li>polynomial division (Newton-based)</li>
<li>polynomial evaluation (Horner) and composition (divide-and-conquer)</li>
<li>product trees, fast multipoint evaluation and interpolation (various algorithms)</li>
<li>power series composition (Horner, Brent-Kung)</li>
<li>added the fmprb_mat module for matrices of balls of real numbers</li>
<li>matrix multiplication</li>
<li>interval-aware LU decomposition, solving, inverse and determinant</li>
<li>many helper functions and small bugfixes</li>
</li>
</ul>
</li>
<li>2012-09-14 - version 0.1</li>
<li>2012-08-05 - began simplified rewrite</li>
<li>2012-04-05 - experimental ball and polynomial code</li>
</ul>

22
doc/introduction.txt
View file

@ -1,22 +0,0 @@
<p>Arb is an experimental C library implementing arbitrary-precision floating-point ball arithmetic, written by Fredrik Johansson &lt;<a href="mailto:fredrik.johansson@gmail.com">fredrik.johansson@gmail.com</a>>. The git repository is <a href="https://github.com/fredrik-johansson/arb/">https://github.com/fredrik-johansson/arb/</a></p>
<p>Ball arithmetic, also known as mid-rad interval arithmetic, is an extension of floating-point arithmetic in which an error bound is attached to each variable. This allows doing rigorous computations over the real numbers, while avoiding the overhead of traditional (inf-sup) interval arithmetic at high precision, and eliminating much of the need for time-consuming and bug-prone manual error analysis.</p>
<p>At the moment, Arb contains:</p>
<ul>
<li>A module (fmpr) for correctly rounded arbitrary-precision floating-point arithmetic. Arb numbers have a few special features, such as arbitrary-size exponents (useful for combinatorics and asymptotics) and dynamic allocation (facilitating implementation of hybrid integer/floating-point and mixed-precision algorithms).</li>
<li>A module (fmprb) for real ball arithmetic, where a ball is implemented as a pair of fmpr numbers.</li>
<li>Functions for fast high-precision computation of some mathematical constants, based on ball arithmetic.</li>
<li>A module (fmprb_poly) for polynomials or power series over the real numbers, implemented using balls as coefficients, with fast polynomial multiplication.</li>
<li>A rudimentary module (fmprb_mat) for matrices over the real numbers, implemented using balls as coefficients.</li>
</ul>
<p>Planned features include: transcendental functions, more extensive polynomial and matrix functionality, and complex balls (and polynomials and matrices thereof).</p>
<p>Arb uses <a href="http://mpir.org/">MPIR</a> and <a href="http://flintlib.org/">FLINT</a> for the underlying integer arithmetic. The code conventions borrow from FLINT, and the project might get merged back into FLINT when the code stabilizes in the future. Arb also uses <a href="http://mpfr.org">MPFR</a> for some fallback code and for testing purposes.</p>
<p>The current version of Arb implements most of its floating-point arithmetic naively using high-level FLINT types. The speed at low precision is far from optimal, and the memory management can sometimes be wasteful. The internals will be rewritten in the future to fix the inefficiencies, which eventually should make Arb ball arithmetic about as fast as mpz or mpfr arithmetic at any precision.</p>
<p><b>Warning</b>: as this is an early version, any part of the interface is subject to change! Also be aware that there are known and unknown bugs.</p>

41
doc/setup.txt
View file

@ -1,41 +0,0 @@
To compile, test and install Arb, do:
<pre>
./configure &lt;options>
make
make check
make install
</pre>
<p>
Arb requires recent versions of MPIR, MPFR and FLINT.
Currently a source checkout of FLINT from
<a href="https://github.com/fredrik-johansson/flint2">https://github.com/fredrik-johansson/flint2</a>
is required.
If MPIR, MPFR or FLINT is installed in some other location than
the default path /usr/local, pass the
flag --with-mpir=... --with-mpfr=... or --with-flint=... with
the correct path to configure (type ./configure --help to show
more options).</p>
<p>Here is a simple sample program to get started using Arb:</p>
<pre>
#include "fmprb.h"
int main()
{
fmprb_t x;
fmprb_init(x);
fmprb_const_pi(x, 50 * 3.33);
fmprb_printd(x, 50); printf("\n");
fmprb_clear(x);
}
</pre>
The output should be something like the following:
<pre>
3.1415926535897932384626433832795028841971693993751 +/- 4.2764e-50
</pre>

246
doc/source/conf.py Normal file
View file

@ -0,0 +1,246 @@
# -*- coding: utf-8 -*-
#
# Arb documentation build configuration file, created by
# sphinx-quickstart on Thu Oct 11 09:33:44 2012.
#
# This file is execfile()d with the current directory set to its containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
import sys, os
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#sys.path.insert(0, os.path.abspath('.'))
# -- General configuration -----------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = ['sphinx.ext.mathjax']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix of source filenames.
source_suffix = '.rst'
# The encoding of source files.
#source_encoding = 'utf-8-sig'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = u'Arb'
copyright = u'2012, Fredrik Johansson'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = '0.3-git'
# The full version, including alpha/beta/rc tags.
release = '0.3-git'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#language = None
# There are two options for replacing |today|: either, you set today to some
# non-false value, then it is used:
#today = ''
# Else, today_fmt is used as the format for a strftime call.
#today_fmt = '%B %d, %Y'
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = []
# The reST default role (used for this markup: `text`) to use for all documents.
default_role = 'math'
latex_preamble = '\usepackage{amsmath,amssymb}\n'
primary_domain = 'c'
# If true, '()' will be appended to :func: etc. cross-reference text.
#add_function_parentheses = True
# If true, the current module name will be prepended to all description
# unit titles (such as .. function::).
#add_module_names = True
# If true, sectionauthor and moduleauthor directives will be shown in the
# output. They are ignored by default.
#show_authors = False
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
html_theme = 'default'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#html_theme_options = {}
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = []
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".
#html_title = None
# A shorter title for the navigation bar. Default is the same as html_title.
#html_short_title = None
# The name of an image file (relative to this directory) to place at the top
# of the sidebar.
#html_logo = None
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
#html_favicon = None
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
#html_last_updated_fmt = '%b %d, %Y'
# If true, SmartyPants will be used to convert quotes and dashes to
# typographically correct entities.
#html_use_smartypants = True
# Custom sidebar templates, maps document names to template names.
#html_sidebars = {}
# Additional templates that should be rendered to pages, maps page names to
# template names.
#html_additional_pages = {}
# If false, no module index is generated.
#html_domain_indices = True
# If false, no index is generated.
#html_use_index = True
# If true, the index is split into individual pages for each letter.
#html_split_index = False
# If true, links to the reST sources are added to the pages.
#html_show_sourcelink = True
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
#html_show_sphinx = True
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
#html_show_copyright = True
# If true, an OpenSearch description file will be output, and all pages will
# contain a <link> tag referring to it. The value of this option must be the
# base URL from which the finished HTML is served.
#html_use_opensearch = ''
# This is the file name suffix for HTML files (e.g. ".xhtml").
#html_file_suffix = None
# Output file base name for HTML help builder.
htmlhelp_basename = 'Arbdoc'
# -- Options for LaTeX output --------------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents = [
('index', 'Arb.tex', u'Arb Documentation',
u'Fredrik Johansson', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
# the title page.
#latex_logo = None
# For "manual" documents, if this is true, then toplevel headings are parts,
# not chapters.
#latex_use_parts = False
# If true, show page references after internal links.
#latex_show_pagerefs = False
# If true, show URL addresses after external links.
#latex_show_urls = False
# Documents to append as an appendix to all manuals.
#latex_appendices = []
# If false, no module index is generated.
#latex_domain_indices = True
# -- Options for manual page output --------------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
('index', 'arb', u'Arb Documentation',
[u'Fredrik Johansson'], 1)
]
# If true, show URL addresses after external links.
#man_show_urls = False
# -- Options for Texinfo output ------------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
('index', 'Arb', u'Arb Documentation',
u'Fredrik Johansson', 'Arb', 'One line description of project.',
'Miscellaneous'),
]
# Documents to append as an appendix to all manuals.
#texinfo_appendices = []
# If false, no module index is generated.
#texinfo_domain_indices = True
# How to display URL addresses: 'footnote', 'no', or 'inline'.
#texinfo_show_urls = 'footnote'

67
doc/source/credits.rst Normal file
View file

@ -0,0 +1,67 @@
Credits and references
===============================================================================
Arb is licensed GNU General Public License version 2, or any later version.
Fredrik's work on Arb is supported by Austrian Science Fund FWF Grant Y464-N18
(Fast Computer Algebra for Special Functions).
Arb includes code by Bill Hart and
Sebastian Pancratz taken from FLINT (also licensed GPL 2.0+).
Software
-------------------------------------------------------------------------------
The following software has been helpful in the development of Arb.
* GMP (Torbjörn Granlund and others), http://gmplib.org
* MPIR (Brian Gladman, Jason Moxham, William Hart and others), http://mpir.org
* MPFR (Guillaume Hanrot, Vincent Lefèvre, Patrick Pélissier, Philippe Théveny, Paul Zimmermann and others), http://mpfr.org
* FLINT (William Hart, Sebastian Pancratz, Andy Novocin, Fredrik Johansson, David Harvey and others), http://flintlib.org
* Sage (William Stein and others), http://sagemath.org
* SymPy (Ondřej Čertík, Aaron Meurer and others), http://sympy.org
* mpmath (Fredrik Johansson and others), http://mpmath.org
* Mathematica (Wolfram Research), http://www.wolfram.com/mathematica
* HolonomicFunctions (Christoph Koutschan), http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/
* Sphinx (George Brandl and others), http://sphinx.pocoo.org
Bibliography
-------------------------------------------------------------------------------
.. [BBC2000] J.Borwein, D. M. Bradley and R. E. Crandall, "Computational strategies for the Riemann zeta function", Journal of Computational and Applied Mathematics 121 (2000) 247-296
.. [BZ1992]_ J.Borwein and I.Zucker, "Fast evaluation of the gamma function for small rational fractions using complete elliptic integrals of the first kind", IMA Journal of Numerical Analysis 12 (1992) 519-526
.. [Bor1987]_ P.Borwein, "Reduced complexity evaluation of hypergeometric functions", Journal of Approximation Theory 50:3 (1987)
.. [Bor2000] P.Borwein, "An Efficient Algorithm for the Riemann Zeta Function", Constructive experimental and nonlinear analysis, CMS Conference Proc. 27 (2000) 29-34. http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf
.. [BM1980] R.Brent and E.M. McMillan, "Some new algorithms for high-precision computation of Euler's constant", Mathematics of Computation 34 (1980) 305-312.
.. [Bre2010] R.Brent, "Ramanujan and Euler's Constant", http://wwwmaths.anu.edu.au/~brent/pd/Euler_CARMA_10.pdf
.. [BZ2011] R.Brent and P.Zimmermann, *Modern Computer Arithmetic*, Cambridge University Press (2011), http://www.loria.fr/~zimmerma/mca/pub226.html
.. [CP2005] R.Crandall and C.Pomerance, *Prime Numbers: A Computational Perspective*, second edition, Springer (2005).
.. [Fil1992] S.Fillebrown, "Faster Computation of Bernoulli Numbers", Journal of Algorithms 13 (1992) 431-445
.. [GG2003] J.von zur Gathen and J. Gerhard, *Modern Computer Algebra*, second edition, Cambridge University Press (2003)
.. [GS2003] X.Gourdon and P. Sebah, "Numerical evaluation of the Riemann Zeta-function" (2003), http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf
.. [HZ2004] G.Hanrot and P. Zimmermann, "Newton Iteration Revisited" (2004), http://www.loria.fr/~zimmerma/papers/fastnewton.ps.gz
.. [Hoe2009] J.van der Hoeven, "Ball arithmetic", Technical Report, HAL 00432152 (2009). http://www.texmacs.org/joris/ball/ball-abs.html
.. [Kar1998] E.A.Karatsuba, "Fast evaluation of the Hurwitz zeta function and Dirichlet L-series", Problems of Information Transmission 34:4 (1998), 342-353. http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=ppi&paperid=425&option_lang=eng
.. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
.. [NIST2012], National Institute of Standards and Technology, *Digital Library of Mathematical Functions* (2012), http://dlmf.nist.gov/
.. [PS1973] M.S. Paterson and L. J. Stockmeyer, "On the number of nonscalar multiplications necessary to evaluate polynomials", SIAM J. Comput (1973)
.. [Tak2000] D.Takahashi, "A fast algorithm for computing large Fibonacci numbers", Information Processing Letters 75 (2000) 243-246. http://www.ii.uni.wroc.pl/~lorys/IPL/article75-6-1.pdf

461
doc/source/fmpr.rst Normal file
View file

@ -0,0 +1,461 @@
fmpr.h -- binary floating-point numbers
===============================================================================
A variable of type *fmpr_t* holds an arbitrary-precision binary
floating-point number, i.e. a rational number of the form
`x \times 2^y` where `x, y \in \mathbb{Z}` and `x` is odd;
or one of the special values zero, plus infinity, minus infinity,
or NaN (not-a-number).
The component `x` is called the *mantissa*, and `y` is called the
*exponent*. Note that this is just one among many possible
conventions: the mantissa (alternatively *significand*) is
sometimes viewed as a fraction in the interval `[1/2, 1)`, with the
exponent pointing to the position above the top bit rather than the
position of the bottom bit, and with a separate sign.
The conventions for special values largely follow those of the
IEEE floating-point standard. At the moment, there is no support
for negative zero, unsigned infinity, or a NaN with a payload, though
some these might be added in the future.
An *fmpr* number is exact and has no inherent "accuracy". We
use the term *precision* to denote either the target precision of
an operation, or the bit size of a mantissa (which in general is
unrelated to the "accuracy" of the number: for example, the
floating-point value 1 has a precision of 1 bit in this sense and is
simultaneously an infinitely accurate approximation of the
integer 1 and a 2-bit accurate approximation of
`\sqrt 2 = 1.011010100\ldots_2`).
Except where otherwise noted, the output of an operation is the
floating-point number obtained by taking the inputs as exact numbers,
in principle carrying out the operation exactly, and rounding the
resulting real number to the nearest representable floating-point
number whose mantissa has at most the specified number of bits, in
the specified direction of rounding. Some operations are always
or optionally done exactly.
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: fmpr_struct
An *fmpr_struct* holds a mantissa and an exponent.
If the mantissa and exponent are sufficiently small, their values are
stored as immediate values in the *fmpr_struct*; large values are
represented by pointers to heap-allocated arbitrary-precision integers.
Currently, both the mantissa and exponent are implemented using
the FLINT *fmpz* type. Special values are currently encoded
by the mantissa being set to zero.
.. type:: fmpr_t
An *fmpr_t* is defined as an array of length one of type
*fmpr_struct*, permitting an *fmpr_t* to be passed by
reference.
.. type:: fmpr_rnd_t
Specifies the rounding mode for the result of an approximate operation.
.. macro:: FMPR_RND_NEAREST
Specifies that the result of an operation should be rounded to the
nearest representable number, rounding to an odd mantissa if there is a tie
between two values. Note: the code for this rounding mode is currently
not implemented.
.. macro:: FMPR_RND_DOWN
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards zero.
.. macro:: FMPR_RND_UP
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction away from zero.
.. macro:: FMPR_RND_FLOOR
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards minus infinity.
.. macro:: FMPR_RND_CEIL
Specifies that the result of an operation should be rounded to the
nearest representable number in the direction towards plus infinity.
.. macro:: FMPR_PREC_EXACT
If passed as the precision parameter to a function, indicates that no
rounding is to be performed. This must only be used when it is known
that the result of the operation can be represented exactly and fits
in memory (the typical use case is working with values small integers).
Note that, for example, adding two numbers whose exponents are far
apart can easily produce an exact result that is far too large to
store in memory.
Memory management
-------------------------------------------------------------------------------
.. function:: void fmpr_init(fmpr_t x)
Initializes the variable *x* for use. Its value is set to zero.
.. function:: void fmpr_clear(fmpr_t x)
Clears the variable *x*, freeing or recycling its allocated memory.
Special values
-------------------------------------------------------------------------------
.. function:: void fmpr_zero(fmpr_t x)
.. function:: void fmpr_one(fmpr_t x)
.. function:: void fmpr_pos_inf(fmpr_t x)
.. function:: void fmpr_neg_inf(fmpr_t x)
.. function:: void fmpr_nan(fmpr_t x)
Sets *x* respectively to 0, 1, `+\infty`, `-\infty`, NaN.
.. function:: int fmpr_is_zero(const fmpr_t x)
.. function:: int fmpr_is_one(const fmpr_t x)
.. function:: int fmpr_is_pos_inf(const fmpr_t x)
.. function:: int fmpr_is_neg_inf(const fmpr_t x)
.. function:: int fmpr_is_nan(const fmpr_t x)
Returns nonzero iff *x* respectively equals
0, 1, `+\infty`, `-\infty`, NaN.
.. function:: int fmpr_is_inf(const fmpr_t x)
Returns nonzero iff *x* equals either `+\infty` or `-\infty`.
.. function:: int fmpr_is_normal(const fmpr_t x)
Returns nonzero iff *x* is a finite, nonzero floating-point value, i.e.
not one of the special values 0, `+\infty`, `-\infty`, NaN.
.. function:: int fmpr_is_special(const fmpr_t x)
Returns nonzero iff *x* is one of the special values
0, `+\infty`, `-\infty`, NaN, i.e. not a finite, nonzero
floating-point value.
Assignment and rounding
-------------------------------------------------------------------------------
.. function:: long _fmpr_normalise(fmpz_t man, fmpz_t exp, long prec, fmpr_rnd_t rnd)
Rounds the mantissa and exponent in-place.
.. function:: void fmpr_set(fmpr_t y, const fmpr_t x)
Sets *y* to a copy of *x*.
.. function:: long fmpr_set_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *y* to a copy of *x* rounded in the direction specified by rnd to the
number of bits specified by prec.
.. function:: void fmpr_set_error_result(fmpr_t err, const fmpr_t result, long rret)
Given the return value *rret* and output variable *result* from a
function performing a rounding (e.g. *fmpr_set_round* or *fmpr_add*), sets
*err* to a bound for the absolute error.
.. function:: void fmpr_add_error_result(fmpr_t err, const fmpr_t err_in, const fmpr_t result, long rret, long prec, fmpr_rnd_t rnd)
Like *fmpr_set_error_result*, but adds *err_in* to the error.
Comparisons
-------------------------------------------------------------------------------
.. function:: int fmpr_equal(const fmpr_t x, const fmpr_t y)
Returns nonzero iff *x* and *y* are exactly equal. This function does
not treat NaN specially, i.e. NaN compares as equal to itself.
.. function:: int fmpr_cmp(const fmpr_t x, const fmpr_t y)
Returns negative, zero, or positive, depending on whether *x* is
respectively smaller, equal, or greater compared to *y*.
Comparison with NaN is undefined.
.. function:: int fmpr_cmpabs(const fmpr_t x, const fmpr_t y)
Compares the absolute values of *x* and *y*.
.. function:: int fmpr_sgn(const fmpr_t x)
Returns `-1`, `0` or `+1` according to the sign of *x*. The sign
of NaN is undefined.
Random number generation
-------------------------------------------------------------------------------
.. function:: void fmpr_randtest(fmpr_t x, flint_rand_t state, long bits, long mag_bits)
Generates a finite random number whose mantissa has precision at most
*bits* and whose exponent has at most *mag_bits* bits. The
values are distributed non-uniformly: special bit patterns are generated
with high probability in order to allow the test code to exercise corner
cases.
.. function:: void fmpr_randtest_not_zero(fmpr_t x, flint_rand_t state, long bits, long mag_bits)
Identical to *fmpr_randtest*, except that zero is never produced
as an output.
.. function:: void fmpr_randtest_special(fmpr_t x, flint_rand_t state, long bits, long mag_bits)
Indentical to *fmpr_randtest*, except that the output occasionally
is set to an infinity or NaN.
Conversions
-------------------------------------------------------------------------------
.. function:: int fmpr_get_mpfr(mpfr_t x, const fmpr_t y, mpfr_rnd_t rnd)
Sets the MPFR variable *x* to the value of *y*. If the
precision of *x* is too small to allow *y* to be represented
exactly, it is rounded in the specified MPFR rounding mode.
The return value indicates the direction of rounding,
following the standard convention of the MPFR library.
.. function:: void fmpr_set_mpfr(fmpr_t x, const mpfr_t y)
Sets *x* to the exact value of the MPFR variable *y*.
.. function:: void fmpr_set_ui(fmpr_t x, ulong c)
.. function:: void fmpr_set_si(fmpr_t x, long c)
.. function:: void fmpr_set_fmpz(fmpr_t x, const fmpz_t c)
Sets *x* exactly to the integer *c*.
.. function:: void fmpr_get_fmpq(fmpq_t y, const fmpr_t x);
Sets *y* to the exact value of *x*. The result is undefined
if *x* is not a finite fraction.
.. function:: long fmpr_set_fmpq(fmpr_t x, const fmpq_t y, long prec, fmpr_rnd_t rnd);
Sets *x* to the value of *y*, rounded according to *prec* and *rnd*.
.. function:: void fmpr_set_fmpz_2exp(fmpr_t x, const fmpz_t man, const fmpz_t exp)
.. function:: void fmpr_set_si_2exp_si(fmpr_t x, long man, long exp)
.. function:: void fmpr_set_ui_2exp_si(fmpr_t x, ulong man, long exp)
Sets *x* to `\mathrm{man} \times 2^{\mathrm{exp}}`.
.. function:: void fmpr_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const fmpr_t x);
Sets *man* and *exp* to the unique integers such that
`x = \mathrm{man} \times 2^{\mathrm{exp}}` and *man* is odd,
provided that *x* is a nonzero finite fraction.
If *x* is zero, both *man* and *exp* are set to zero. If *x* is
infinite or NaN, the result is undefined.
.. function:: int fmpr_get_fmpz_fixed_fmpz(fmpz_t y, const fmpr_t x, const fmpz_t e)
.. function:: int fmpr_get_fmpz_fixed_si(fmpz_t y, const fmpr_t x, long e)
Converts *x* to a mantissa with predetermined exponent, i.e. computes
an integer *y* such that `y \times 2^e \approx x`, truncating if necessary.
Returns 0 if exact and 1 if truncation occurred.
Input and output
-------------------------------------------------------------------------------
.. function:: void fmpr_print(const fmpr_t x)
Prints the mantissa and exponent of *x* as integers, precisely showing
the internal representation.
.. function:: void fmpr_printd(const fmpr_t x, long digits)
Prints *x* as a decimal floating-point number, rounding to the specified
number of digits. This function is currently implemented using MPFR,
and does not support large exponents.
Arithmetic
-------------------------------------------------------------------------------
.. function:: void fmpr_neg(fmpr_t y, const fmpr_t x)
Sets *y* to the negation of *x*.
.. function:: long fmpr_neg_round(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *y* to the negation of *x*, rounding the result.
.. function:: void fmpr_abs(fmpr_t y, const fmpr_t x)
Sets *y* to the absolute value of *x*.
.. function:: long fmpr_add(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_add_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_add_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_add_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = x + y`, rounded according to *prec* and *rnd*. The precision
can be *FMPR_PREC_EXACT* to perform an exact addition, provided that the
result fits in memory.
.. function:: long _fmpr_add_eps(fmpr_t z, const fmpr_t x, int sign, long prec, fmpr_rnd_t rnd)
Sets *z* to the value that results by adding an infinitesimal quantity
of the given sign to *x*, and rounding. The result is undefined
if *x* is zero.
.. function:: long fmpr_sub(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_sub_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_sub_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_sub_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = x - y`, rounded according to *prec* and *rnd*. The precision
can be *FMPR_PREC_EXACT* to perform an exact addition, provided that the
result fits in memory.
.. function:: long fmpr_mul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_mul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_mul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_mul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = x \times y`, rounded according to prec and rnd. The precision
can be *FMPR_PREC_EXACT* to perform an exact multiplication, provided that the
result fits in memory.
.. function:: void fmpr_mul_2exp_si(fmpr_t y, const fmpr_t x, long e)
.. function:: void fmpr_mul_2exp_fmpz(fmpr_t y, const fmpr_t x, const fmpz_t e)
Sets *y* to *x* multiplied by `2^e` without rounding.
.. function:: long fmpr_div(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_div_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_ui_div(fmpr_t z, ulong x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_div_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_si_div(fmpr_t z, long x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_div_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_fmpz_div(fmpr_t z, const fmpz_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_fmpz_div_fmpz(fmpr_t z, const fmpz_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = x / y`, rounded according to *prec* and *rnd*. If *y* is zero,
*z* is set to NaN.
.. function:: long fmpr_addmul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_addmul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_addmul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_addmul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = z + x \times y`, rounded according to *prec* and *rnd*. The
intermediate multiplication is always performed without roundoff. The
precision can be *FMPR_PREC_EXACT* to perform an exact addition, provided
that the result fits in memory.
.. function:: long fmpr_submul(fmpr_t z, const fmpr_t x, const fmpr_t y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_submul_ui(fmpr_t z, const fmpr_t x, ulong y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_submul_si(fmpr_t z, const fmpr_t x, long y, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_submul_fmpz(fmpr_t z, const fmpr_t x, const fmpz_t y, long prec, fmpr_rnd_t rnd)
Sets `z = z - x \times y`, rounded according to *prec* and *rnd*. The
intermediate multiplication is always performed without roundoff. The
precision can be *FMPR_PREC_EXACT* to perform an exact subtraction, provided
that the result fits in memory.
.. function:: long fmpr_sqrt(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_sqrt_ui(fmpr_t z, ulong x, long prec, fmpr_rnd_t rnd)
.. function:: long fmpr_sqrt_fmpz(fmpr_t z, const fmpz_t x, long prec, fmpr_rnd_t rnd)
Sets *z* to the square root of *x*, rounded according to *prec* and *rnd*.
The result is NaN if *x* is negative.
.. function:: void fmpr_pow_sloppy_fmpz(fmpr_t y, const fmpr_t b, const fmpz_t e, long prec, fmpr_rnd_t rnd)
.. function:: void fmpr_pow_sloppy_ui(fmpr_t y, const fmpr_t b, ulong e, long prec, fmpr_rnd_t rnd)
.. function:: void fmpr_pow_sloppy_si(fmpr_t y, const fmpr_t b, long e, long prec, fmpr_rnd_t rnd)
Sets `y = b^e`, computed using without guaranteeing correct (optimal)
rounding, but guaranteeing that the result is a correct upper or lower
bound if the rounding is directional. Currently requires `b \ge 0`.
Special functions
-------------------------------------------------------------------------------
.. function:: long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *z* to `\log(x)`, rounded according to *prec* and *rnd*.
The result is NaN if *x* is negative.
This function is currently implemented using MPFR and does not
support large exponents.
.. function:: long fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *z* to `\log(1+x)`, rounded according to *prec* and *rnd*.
This function
computes an accurate value when *x* is small.
The result is NaN if `1+x` is negative.
This function is currently implemented using MPFR and does not
support large exponents.
.. function:: long fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *z* to `\exp(x)`, rounded according to *prec* and *rnd*.
This function is currently implemented using MPFR and does not
support large exponents.
.. function:: long fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
Sets *z* to `\exp(x)-1`, rounded according to *prec* and *rnd*.
This function computes an accurate value when *x* is small.
This function is currently implemented using MPFR and does not
support large exponents.

662
doc/source/fmprb.rst Normal file
View file

@ -0,0 +1,662 @@
fmprb.h -- real numbers represented as floating-point balls
===============================================================================
An *fmprb_t* represents a ball over the real numbers,
that is, an interval `[m-r, m+r]` where the midpoint `m` and the
radius `r` are (extended) real numbers and `r` is nonnegative.
The result of an (approximate) operation done on *fmprb_t* variables
is a ball which contains the result of the (mathematically exact) operation
applied to any choice of points in the input balls.
In general, the output ball is not the smallest possible.
The precision parameter passed to each function roughly indicates the
precision to which calculations on the midpoint are carried out
(operations on the radius are always done using a fixed, small
precision.)
For arithmetic operations, the precision parameter currently simply
specifies the precision of the corresponding *fmpr* operation.
In the future, the arithmetic might be made faster by incorporating
sloppy rounding (typically equivalent to a loss of 1-2 bits of effective
working precision) when the result is known to be inexact (while still
propagating errors rigorously, of course).
Arithmetic operations done on exact input with exactly
representable output are always guaranteed to produce exact output.
For more complex operations, the precision parameter indicates a minimum
working precision (algorithms might allocate extra internal precision to
attempt to produce an output accurate to the requested number of bits,
especially when the required precision can be estimated easily, but this
is not generally required).
If the precision is increased and the inputs either are exact or are
computed with increased accuracy as well, the output should
converge proportionally, absent any bugs.
The general intended strategy for using ball arithmetic is to add a few
guard bits, and then repeat the calculation as necessary with an
exponentially increasing number of guard bits (Ziv's strategy) until the
result is exact
enough for one's purposes (typically the first attempt will be successful).
There are some caveats: in general, ball arithmetic only makes
sense for (Lipschitz) continuous functions, and
trying to approximate functions close to singularities might result in
slow convergence, or failure to converge.
**Warning**: some methods for transcendental functions and constants
currently perform the error propagation in a non-rigorous way, due to the
implementation being incomplete (in some cases, a
rigorous error bound for the algorithm or function might not
be known at all).
This should be indicated in the documentation for each function.
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: fmprb_struct
.. type:: fmprb_t
An *fmprb_struct* consists of a pair of *fmpr_struct*:s.
An *fmprb_t* is defined as an array of length one of type
*fmprb_struct*, permitting an *fmprb_t* to be passed by
reference.
.. macro:: FMPRB_RAD_PREC
The precision used for operations on the radius. This is small
enough to fit in a single word, currently 30 bits.
.. macro:: fmprb_midref(x)
Macro returning a pointer to the midpoint of *x* as an *fmpr_t*.
.. macro:: fmprb_radref(x)
Macro returning a pointer to the radius of *x* as an *fmpr_t*.
Memory management
-------------------------------------------------------------------------------
.. function:: void fmprb_init(fmprb_t x)
Initializes the variable *x* for use. Its midpoint and radius are both
set to zero.
.. function:: void fmprb_clear(fmprb_t x)
Clears the variable *x*, freeing or recycling its allocated memory.
.. function:: fmprb_struct * _fmprb_vec_init(long n)
Returns a pointer to an array of *n* initialized *fmprb_struct*:s.
.. function:: void _fmprb_vec_clear(fmprb_struct * v, long n)
Clears an array of *n* initialized *fmprb_struct*:s.
Basic manipulation
-------------------------------------------------------------------------------
.. function:: int fmprb_is_exact(const fmprb_t x)
Returns nonzero iff the radius of *x* is zero.
.. function:: int fmprb_equal(const fmprb_t x, const fmprb_t y)
Returns nonzero iff *x* and *y* are equal as balls, i.e. have both the
same midpoint and radius.
.. function:: void fmprb_zero(fmprb_t x)
Sets *x* to zero.
.. function:: int fmprb_is_zero(const fmprb_t x)
Returns nonzero iff the midpoint and radius of *x* are both zero.
.. function:: void fmprb_set(fmprb_t y, const fmprb_t x)
Sets *y* to a copy of *x*.
.. function:: void fmprb_set_round(fmprb_t y, const fmprb_t x, long prec)
Sets *y* to a copy of *x*, rounded to *prec* bits.
.. function:: void fmprb_neg(fmprb_t y, const fmprb_t x)
Sets *y* to the negation of *x*.
.. function:: void fmprb_abs(fmprb_t y, const fmprb_t x)
Sets *y* to the absolute value of *x*. No attempt is made to improve the
interval represented by *x* if it contains zero.
.. function:: void fmprb_set_fmpr(fmprb_t y, const fmpr_t x)
.. function:: void fmprb_set_si(fmprb_t y, long x)
.. function:: void fmprb_set_ui(fmprb_t y, ulong x)
.. function:: void fmprb_set_fmpz(fmprb_t y, const fmpz_t x)
Sets *y* exactly to *x*.
.. function:: void fmprb_set_fmpq(fmprb_t y, const fmpq_t x, long prec)
Sets *y* to the rational number *x*, rounded to *prec* bits.
.. function:: void fmprb_one(fmprb_t x)
Sets *x* to the exact integer 1.
.. function:: int fmprb_is_one(const fmprb_t x)
Returns nonzero iff *x* is exactly 1.
Input and output
-------------------------------------------------------------------------------
.. function:: void fmprb_print(const fmprb_t x)
Prints the internal representation of *x*.
.. function:: void fmprb_printd(const fmprb_t x, long digits)
Prints *x* in decimal. The printed value of the radius is not adjusted
to compensate for the fact that the binary-to-decimal conversion
of both the midpoint and the radius introduces additional error.
Random number generation
-------------------------------------------------------------------------------
.. function:: void fmprb_randtest(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
Generates a random ball. The midpoint and radius will both be finite.
.. function:: void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, long bits)
Sets *q* to a random rational number from the interval represented by *x*.
A denominator is chosen by multiplying the binary denominator of *x*
by a random integer up to *bits* bits.
The outcome is undefined if the midpoint or radius of *x* is non-finite,
or if the exponent of the midpoint or radius is so large or small
that representing the endpoints as exact rational numbers would
cause overflows.
Precision and comparisons
-------------------------------------------------------------------------------
.. function:: void fmprb_add_error_fmpr(fmprb_t x, const fmpr_t err)
Adds *err*, which is assumed to be nonnegative, to the radius of *x*.
.. function:: void fmprb_add_error_2exp_si(fmprb_t x, long e)
Adds `2^e` to the radius of *x*.
.. function:: void fmprb_add_error(fmprb_t x, const fmprb_t err)
Adds the supremum of *err*, which is assumed to be nonnegative, to the
radius of *x*.
.. function:: int fmprb_contains_fmpr(const fmprb_t x, const fmpr_t y)
.. function:: int fmprb_contains_fmpq(const fmprb_t x, const fmpq_t y)
.. function:: int fmprb_contains_fmpz(const fmprb_t x, const fmpz_t y)
.. function:: int fmprb_contains_mpfr(const fmprb_t x, const mpfr_t y)
.. function:: int fmprb_contains_zero(const fmprb_t x)
Returns nonzero iff the given number is contained in the interval
represented by *x*.
.. function:: void fmprb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const fmprb_t x)
Computes the exact interval represented by *x*, in the form of an integer
interval multiplied by a power of two, i.e. `x = [a, b] \times 2^{\mathrm{exp}}`.
The outcome is undefined if the midpoint or radius of *x* is non-finite,
or if the difference in magnitude between the midpoint and radius
is so large that representing the endpoints exactly would cause overflows.
.. function:: long fmprb_rel_error_bits(const fmprb_t x)
Returns the effective relative error of *x* measured in bits, defined as
the difference between the position of the top bit in the radius
and the top bit in the midpoint, plus one.
The result is clamped between plus/minus *FMPR_PREC_EXACT*.
.. function:: long fmprb_rel_accuracy_bits(const fmprb_t x)
Returns the effective relative accuracy of *x* measured in bits,
equal to the negative of the return value from *fmprb_rel_error_bits*.
Arithmetic
-------------------------------------------------------------------------------
.. function:: void fmprb_add(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_add_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_add_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_add_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
.. function:: void fmprb_add_fmpr(fmprb_t z, const fmprb_t x, const fmpr_t y, long prec)
Sets `z = x + y`, rounded to *prec* bits. The precision can be
*FMPR_PREC_EXACT* provided that the result fits in memory.
.. function:: void fmprb_sub(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_sub_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_sub_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_sub_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets `z = x - y`, rounded to *prec* bits. The precision can be
*FMPR_PREC_EXACT* provided that the result fits in memory.
.. function:: void fmprb_mul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_mul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_mul_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_mul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets `z = x \times y`, rounded to *prec* bits. The precision can be
*FMPR_PREC_EXACT* provided that the result fits in memory.
.. function:: void fmprb_mul_2exp_si(fmprb_t y, const fmprb_t x, long e)
Sets *y* to *x* multiplied by `2^e`.
.. function:: void fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_div_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_div_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_div_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
.. function:: void fmprb_fmpz_div_fmpz(fmprb_t y, const fmpz_t num, const fmpz_t den, long prec)
.. function:: void fmprb_ui_div(fmprb_t z, ulong x, const fmprb_t y, long prec);
Sets `z = x / y`, rounded to *prec* bits. If *y* contains zero, *z* is
set to `0 \pm \infty`. Otherwise, error propagation uses the rule
.. math ::
\left| \frac{x}{y} - \frac{x+\xi_1 a}{y+\xi_2 b} \right| =
\left|\frac{x \xi_2 b - y \xi_1 a}{y (y+\xi_2 b)}\right| \le
\frac{|xb|+|ya|}{|y| (|y|-b)}
where `-1 \le \xi_1, \xi_2 \le 1`, and
where the triangle inequality has been applied to the numerator and
the reverse triangle inequality has been applied to the denominator.
.. function:: void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec);
Sets `y = x / (2^n - 1)`, rounded to *prec* bits.
.. function:: void fmprb_addmul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_addmul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_addmul_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_addmul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets `z = z + x \times y`, rounded to prec bits. The precision can be
*FMPR_PREC_EXACT* provided that the result fits in memory.
.. function:: void fmprb_submul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_submul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)
.. function:: void fmprb_submul_si(fmprb_t z, const fmprb_t x, long y, long prec)
.. function:: void fmprb_submul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, long prec)
Sets `z = z - x \times y`, rounded to *prec* bits. The precision can be
*FMPR_PREC_EXACT* provided that the result fits in memory.
.. function:: void fmprb_sqrt(fmprb_t z, const fmprb_t x, long prec)
.. function:: void fmprb_sqrt_ui(fmprb_t z, ulong x, long prec)
.. function:: void fmprb_sqrt_fmpz(fmprb_t z, const fmpz_t x, long prec)
Sets *z* to the square root of *x*, rounded to *prec* bits.
Error propagation is done using the following rule:
assuming `m > r \ge 0`, the error is largest at `m - r`, and we have
`\sqrt{m} - \sqrt{m-r} \le r / (2 \sqrt{m-r})`.
void fmprb_pow_fmpz(fmprb_t y, const fmprb_t b, const fmpz_t e, long prec)
.. function:: void fmprb_pow_ui(fmprb_t y, const fmprb_t b, ulong e, long prec)
.. function:: void fmprb_ui_pow_ui(fmprb_t y, ulong b, ulong e, long prec)
.. function:: void fmprb_si_pow_ui(fmprb_t y, long b, ulong e, long prec)
Sets `y = b^e` using binary exponentiation. Provided that *b* and *e*
are small enough and the exponent is positive, the exact power can be
computed using *FMPR_PREC_EXACT*.
Special functions
-------------------------------------------------------------------------------
Elementary functions
...............................................................................
.. function:: void fmprb_log(fmprb_t z, const fmprb_t x, long prec)
.. function:: void fmprb_log_ui(fmprb_t z, ulong x, long prec)
.. function:: void fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec)
Sets `z = \log(x)`. Error propagation is done using the following rule:
assuming `m > r \ge 0`, the error is largest at `m - r`, and we have
`\log(m) - \log(m-r) = \log(1 + r/(m-r))`. The last expression is
calculated accurately for small radii via *fmpr_log1p*.
An input containing zero currently raises an exception.
.. function:: void fmprb_exp(fmprb_t z, const fmprb_t x, long prec)
Sets `z = \exp(x)`. Error propagation is done using the following rule:
the error is largest at `m + r`, and we have
`\exp(m+r) - \exp(m) = \exp(m) (\exp(r)-1)`.
The last expression is calculated accurately for small radii
via *fmpr_expm1*.
.. function:: void fmprb_sin(fmprb_t s, const fmprb_t x, long prec)
.. function:: void fmprb_cos(fmprb_t c, const fmprb_t x, long prec)
.. function:: void fmprb_sin_cos(fmprb_t s, fmprb_t c, const fmprb_t x, long prec)
Sets `s = \sin x`, `c = \cos x`. Error propagation uses the rule
`|\sin(m \pm r) - \sin(m)| \le r` (this could be tightened to
`\min(r,2)`).
.. function:: void fmprb_atan(fmprb_t z, const fmprb_t x, long prec)
Sets `z = \tan^{-1} x`. Letting `d = \max(0, |m| - r)`,
the propagated error is bounded by `r / (1 + d^2)`
(this could be tightened).
Factorials and other integer functions
...............................................................................
.. function:: void fmprb_fac_ui(fmprb_t x, ulong n, long prec)
Sets *x* to *n* factorial, computed using binary splitting. Provided that
*n* is small enough, the exact factorial can be computed using
*FMPR_PREC_EXACT*.
.. function:: void fmprb_rfac_ui_bsplit(fmprb_t y, const fmprb_t x, ulong n, long prec)
Sets *x* to the rising factorial `x (x+1) (x+2) \cdots (x+n-1)`,
computed using binary splitting.
The basecase processes eight factors at a time using the formula
.. math ::
x(x+1)\cdots(x+7) = (28 + 98x + 63x^2 + 14x^3 + x^4)^2 - 16 (7+2x)^2,
replacing 7 full-precision multiplications with 4 squarings and
1 multiplication ([CP2005]_, page 316).
Empirically, this is about twice as fast at high precision.
Numerical stability is slightly worse.
.. function:: void fmprb_rfac_ui_multipoint(fmprb_t y, const fmprb_t x, ulong n, long prec)
Sets *x* to the rising factorial `x (x+1) (x+2) \cdots (x+n-1)`,
computed using fast multipoint evaluation. This only requires
`O(n^{1/2+\varepsilon})` multiplications, but has high overhead
and poor numerical stability (adding `O(n)` guard bits to the input
might be necessary to achieve full accuracy). It can be expected to
be faster than the binary splitting algorithm if the input is a
full-precision number, the precision is at least 100000 bits,
and *n* is of the same order of magnitude as (perhaps slightly
smaller than) the number of bits.
.. function:: void fmprb_bin_ui(fmprb_t x, const fmprb_t n, ulong k, long prec)
.. function:: void fmprb_bin_uiui(fmprb_t x, ulong n, ulong k, long prec)
Sets *x* to the binomial coefficient `{n \choose k}`, computed using
binary splitting. Provided that *n* and *k* are small enough, an exact
binomial coefficient can be computed using *FMPR_PREC_EXACT*.
.. function:: void fmprb_fib_fmpz(fmprb_t f, const fmpz_t n, long prec)
.. function:: void fmprb_fib_ui(fmprb_t f, ulong n, long prec)
Sets x to the Fibonacci number `F_n`. Uses the binary squaring
algorithm described in [Tak2000]_.
Provided that *n* is small enough, an exact Fibonacci number can be
computed using *FMPR_PREC_EXACT*.
Constants
...............................................................................
.. function:: void fmprb_const_pi_chudnovsky(fmprb_t x, long prec)
Sets *x* to `\pi`, computed using the Chudnovsky algorithm.
Letting `A = 13591409`, `B = 545140134`, `C = 640320`,
we have `\pi \approx 1 / s_N` where
.. math ::
s_N = 12 \sum_{k=0}^N \frac{(-1)^k (6k)! (A+Bk)}
{(3k)! (k!)^3 C^{3k+3/2}}
The implementation computes an approximation for the
algebraic number `1/s_N` using binary splitting, bounding
the rounding error automatically.
The hypergeometric term ratio is asymptotically
`R = C^3 / (2^6 \times 3^3) \approx 1.5 \times 10^{14}`, and in fact we have
`|\pi - 1/s_N| < 1/R^N` (with a more detailed calculation, the truncation
error could be bounded closer to `1/R^{N+1}`).
.. function:: void fmprb_const_pi(fmprb_t x, long prec)
Sets *x* to `\pi`. The value is cached for repeated use.
.. function:: void fmprb_const_log_sqrt2pi(fmprb_t x, long prec)
Sets *x* to `\log \sqrt{2 \pi}`. The value is cached for repeated use.
.. function:: void fmprb_const_euler_brent_mcmillan(fmprb_t res, long prec)
Sets *x* to Euler's constant `\gamma`, computed using the second
Bessel function formula of Brent and McMillan ([BM1980]_, [MPFR2012]_).
Brent and McMillan conjectured that the error depending
on the internal parameter *n* is of order `O(e^{-8n})`. Brent has
recently proved that this bound is correct, but without determining
an explicit big-O factor [Bre2010]_.
Riemann zeta function
...............................................................................
.. function:: void fmprb_const_zeta3_bsplit(fmprb_t x, long prec)
Sets *x* to Apery's constant `\zeta(3)`, computed by applying binary
splitting to a hypergeometric series.
.. function:: void fmprb_zeta_ui_asymp(fmprb_t z, ulong s, long prec)
Assuming `s \ge 2`, approximates `\zeta(s)` by `1 + 2^{-s}` along with
a correct error bound. We use the following bounds: for `s > b`,
`\zeta(s) - 1 < 2^{-b}`, and generally,
`\zeta(s) - (1 + 2^{-s}) < 2^{2-\lfloor 3 s/2 \rfloor}`.
.. function:: void fmprb_zeta_ui_euler_product(fmprb_t z, ulong s, long prec)
Computes `\zeta(s)` using the Euler product. This is fast only if *s*
is large compared to the precision.
Writing `P(a,b) = \prod_{a \le p \le b} (1 - p^{-s})`, we have
`1/\zeta(s) = P(a,M) P(M+1,\infty)`.
To bound the error caused by truncating
the product at `M`, we write `P(M+1,\infty) = 1 - \epsilon(s,M)`.
Since `0 < P(a,M) \le 1`, the absolute error for `\zeta(s)` is
bounded by `\epsilon(s,M)`.
According to the analysis in [Fil1992]_, it holds for all `s \ge 6` and `M \ge 1`
that `1/P(M+1,\infty) - 1 \le f(s,M) \equiv 2 M^{1-s} / (s/2 - 1)`.
Thus, we have `1/(1-\epsilon(s,M)) - 1 \le f(s,M)`, and expanding
the geometric series allows us to conclude that
`\epsilon(M) \le f(s,M)`.
.. function:: void fmprb_zeta_ui_bernoulli(fmprb_t x, ulong n, long prec)
Computes `\zeta(n)` for even *n* via the corresponding Bernoulli number,
which is generated using FLINT.
.. function:: void fmprb_zeta_ui_vec_borwein(fmprb_struct * z, ulong start, long num, ulong step, long prec)
Evaluates `\zeta(s)` at `\mathrm{num}` consecutive integers *s* beginning
with *start* and proceeding in increments of *step*.
Uses Borwein's formula ([Bor2000]_, [GS2003]_),
implemented to support fast multi-evaluation
(but also works well for a single *s*).
Requires `\mathrm{start} \ge 2`. For efficiency, the largest *s*
should be at most about as
large as *prec*. Arguments approaching *LONG_MAX* will cause
overflows.
One should therefore only use this function for *s* up to about *prec*, and
then switch to the Euler product.
The algorithm for single *s* is basically identical to the one used in MPFR
(see [MPFR2012]_ for a detailed description).
In particular, we evaluate the sum backwards to avoid storing more than one
`d_k` coefficient, and use integer arithmetic throughout since it
is convenient and the terms turn out to be slightly larger than
`2^\mathrm{prec}`.
The only numerical error in the main loop comes from the division by `k^s`,
which adds less than 1 unit of error per term.
For fast multi-evaluation, we repeatedly divide by `k^{\mathrm{step}}`.
Each division reduces the input error and adds at most 1 unit of
additional rounding error, so by induction, the error per term
is always smaller than 2 units.
.. function:: void fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec)
Computes `\zeta(s)` for arbitrary `s \ge 2` using a binary splitting
implementation of Borwein's formula. The algorithm has quasilinear
complexity with respect to the precision.
.. function:: void fmprb_zeta_ui(fmprb_t x, ulong s, long prec)
Computes `\zeta(s)` for nonnegative integer `s \ne 1`, automatically
choosing an appropriate algorithm.
.. function:: void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec)
.. function:: void fmprb_zeta_ui_vec_even(fmprb_struct * x, ulong start, long num, long prec)
.. function:: void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec)
Computes `\zeta(s)` at num consecutive integers (respectively num
even or num odd integers) beginning with `s = \mathrm{start} \ge 2`,
automatically choosing an appropriate algorithm.
Gamma function
...............................................................................
.. function:: void fmprb_gamma_fmpq_karatsuba(fmprb_struct * v, const fmpq_t a, long num, long prec)
Uses Karatsuba's algorithm [Kar1998]_ to compute num coefficients in the
Taylor series of `\Gamma(a+x)` for rational `0 < a \le 1`, i.e.
computes `\Gamma(a), \Gamma'(a) ... \Gamma^{(\mathrm{num}-1)}(a) / (\mathrm{num}-1)!`
This algorithm is most efficient at high precision, for num much smaller
than the number of bits, and with small denominators of *a*.
In particular, with num = 1, this algorithm computes `\Gamma(a)`
efficiently for small rational *a*.
Let `s = \max(2, \mathrm{num}-1)`. With parameters `r` and `n` chosen
such that `r \ge n` and `n \ge 2 s \log 2 s`, Karatsuba shows that
.. math ::
\Gamma^{(j)}(a) = \sum_{k=0}^r \frac{(-1)^k}{k!}
\frac{n^{k+a}}{k+a}
\sum_{m=0}^j (-1)^m \frac{j! \, \log^{j-m} n}{(j-m)! (k+a)^m}
+ \theta_j
where
.. math ::
|\theta_j| \le \frac{5}{3} e^{-n} \log^s n +
\left(\frac{e}{r+2}\right)^{r+2} (1 + n^{r+2} \log^s n).
We choose the parameters `n` and `r` heuristically to be nearly optimal,
and then evaluate the above formula to bound `\theta_j` rigorously.
Karatsuba claims that choosing `r \ge 3n` gives
`|\theta_j| \le 2^{-n-1}`. This is, unfortunately, incorrect.
Setting `r = n \alpha` and expanding the error term around `n = \infty`,
one finds that `\alpha` asymptotically should be
`1/W(1/e) \approx 3.59112147666862` where `W(x)` is the Lambert W-function.
We also optimize the selection of `n` by choosing
`n \approx b \log 2` where `b` is the desired number of bits, rather
than `n \approx b`, and round `n` so that it has a short binary expansion
(this gives smaller numbers in the binary splitting stage).
Finally, if `s` is small, we perform binary splitting to a working
precision of about `2.2` times the target precision rather than exactly.
This factor was tested to give full accuracy up to at least one million
digits when `s \approx 1`. A more careful analysis should
be done here so that a working precision is selected which always is
sufficient and also nearly optimal.
.. function:: void fmprb_gamma_log(fmprb_t y, const fmprb_t x, long prec)
Sets `y = \log \Gamma(x)`, assuming that `x > 0`.
For large `x`, uses Stirling's expansion
.. math ::
\log \Gamma(x) = \left(x-\frac{1}{2}\right)\log x - x +
\frac{\ln {2 \pi}}{2} + \sum_{k=1}^{n-1} t_k + R(n,x)
where
.. math ::
t_k = \frac{B_{2k}}{2k(2k-1)x^{2k-1}}
and `|R(n,x)| < t_n`.
If `x` is too small for the asymptotic expansion to give sufficient
accuracy directly, we translate to `x + r`
using the formula `\log \Gamma(x) = \log \Gamma(x+r) -
\log(x (x+1) (x+2) \cdots (x+r-1))`.
To obtain a remainder smaller than `2^{-b}`, we must choose an `r` such
that `x + r > \beta b`, where `\beta > \log(2) / (2 \pi) \approx 0.11`.
We use a slightly larger factor `\beta \approx 0.2` to more closely
balance `n` and `r`. A much larger `\beta` (e.g. `\beta = 1`) could be
used to reduce the number of Bernoulli numbers that have to be
precomputed, at the expense of slower repeated evaluation.

181
doc/source/fmprb_mat.rst Normal file
View file

@ -0,0 +1,181 @@
fmprb_mat.h -- matrices of real number balls
===============================================================================
An *fmprb_mat_t* represents a dense matrix over the real numbers,
implemented as an array of entries of type *fmprb_struct*.
The dimension (number of rows and columns) of a matrix is fixed at
initialization, and the user must ensure that inputs and outputs to
an operation have compatible dimensions. The number of rows or columns
in a matrix can be zero.
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: fmprb_mat_struct
.. type:: fmprb_mat_t
Contains a pointer to a flat array of the entries (entries), an array of
pointers to the start of each row (rows), and the number of rows (r)
and columns (c).
An *fmprb_mat_t* is defined as an array of length one of type
*fmprb_mat_struct*, permitting an *fmprb_mat_t* to
be passed by reference.
.. macro:: fmprb_mat_entry(mat, i, j)
Macro giving a pointer to the entry at row *i* and column *j*.
.. macro:: fmprb_mat_nrows(mat)
Returns the number of rows of the matrix.
.. macro:: fmprb_mat_ncols(mat)
Returns the number of columns of the matrix.
Memory management
-------------------------------------------------------------------------------
.. function:: void fmprb_mat_init(fmprb_mat_t mat, long r, long c)
Initializes the matrix, setting it to the zero matrix with *r* rows
and *c* columns.
.. function:: void fmprb_mat_clear(fmprb_mat_t mat)
Clears the matrix, deallocating all entries.
Conversions
-------------------------------------------------------------------------------
.. function:: void fmprb_mat_set(fmprb_mat_t dest, const fmprb_mat_t src)
.. function:: void fmprb_mat_set_fmpz_mat(fmprb_mat_t dest, const fmpz_mat_t src)
.. function:: void fmprb_mat_set_fmpq_mat(fmprb_mat_t dest, const fmpq_mat_t src, long prec)
Sets *dest* to *src*. The operands must have identical dimensions.
Input and output
-------------------------------------------------------------------------------
.. function:: void fmprb_mat_printd(const fmprb_mat_t mat, long digits)
Prints each entry in the matrix with the specified number of decimal digits.
Comparisons
-------------------------------------------------------------------------------
.. function:: int fmprb_mat_equal(const fmprb_mat_t mat1, const fmprb_mat_t mat2)
Returns nonzero iff the matrices have the same dimensions
and identical entries.
.. function:: int fmprb_mat_contains_fmpq_mat(const fmprb_mat_t mat1, const fmpq_mat_t mat2)
Returns nonzero iff the matrices have the same dimensions and each entry
in *mat2* is contained in the corresponding entry in *mat1*.
Special matrices
-------------------------------------------------------------------------------
.. function:: void fmprb_mat_zero(fmprb_mat_t mat)
Sets all entries in mat to zero.
.. function:: void fmprb_mat_one(fmprb_mat_t mat)
Sets the entries on the main diagonal to ones,
and all other entries to zero.
Arithmetic
-------------------------------------------------------------------------------
.. function:: void fmprb_mat_neg(fmprb_mat_t dest, const fmprb_mat_t src)
Sets *dest* to the exact negation of *src*. The operands must have
the same dimensions.
.. function:: void fmprb_mat_add(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t mat2, long prec)
Sets res to the sum of *mat1* and *mat2*. The operands must have the same dimensions.
.. function:: void fmprb_mat_sub(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t mat2, long prec)
Sets *res* to the difference of *mat1* and *mat2*. The operands must have
the same dimensions.
.. function:: void fmprb_mat_mul(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t mat2, long prec)
Sets *res* to the matrix product of *mat1* and *mat2*. The operands must have
compatible dimensions for matrix multiplication.
Gaussian elimination and solving
-------------------------------------------------------------------------------
.. function:: int fmprb_mat_lu(long * perm, fmprb_mat_t LU, const fmprb_mat_t A, long prec)
Given an `n \times n` matrix `A`, computes an LU decomposition `PLU = A`
using Gaussian elimination with partial pivoting.
The input and output matrices can be the same, performing the
decomposition in-place.
Entry `i` in the permutation vector perm is set to the row index in
the input matrix corresponding to row `i` in the output matrix.
The algorithm succeeds and returns nonzero if it can find `n` invertible
(i.e. not containing zero) pivot entries. This guarantees that the matrix
is invertible.
The algorithm fails and returns zero, leaving the entries in `P` and `LU`
undefined, if it cannot find `n` invertible pivot elements.
In this case, either the matrix is singular, the input matrix was
computed to insufficient precision, or the LU decomposition was
attempted at insufficient precision.
.. function:: void fmprb_mat_solve_lu_precomp(fmprb_mat_t X, const long * perm, const fmprb_mat_t LU, const fmprb_mat_t B, long prec)
Solves `AX = B` given the precomputed nonsingular LU decomposition `A = PLU`.
The matrices `X` and `B` are allowed to be aliased with each other,
but `X` is not allowed to be aliased with `LU`.
.. function:: int fmprb_mat_solve(fmprb_mat_t X, const fmprb_mat_t A, const fmprb_mat_t B, long prec)
Solves `AX = B` where `A` is a nonsingular `n \times n` matrix
and `X` and `B` are `n \times m` matrices, using LU decomposition.
If `m > 0` and `A` cannot be inverted numerically (indicating either that
`A` is singular or that the precision is insufficient), the values in the
output matrix are left undefined and zero is returned. A nonzero return
value guarantees that `A` is invertible and that the exact solution
matrix is contained in the output.
.. function:: int fmprb_mat_inv(fmprb_mat_t X, const fmprb_mat_t A, long prec)
Sets `X = A^{-1}` where `A` is a square matrix, computed by solving
the system `AX = I`.
If `A` cannot be inverted numerically (indicating either that
`A` is singular or that the precision is insufficient), the values in the
output matrix are left undefined and zero is returned.
A nonzero return value guarantees that the matrix is invertible
and that the exact inverse is contained in the output.
.. function:: void fmprb_mat_det(fmprb_t det, const fmprb_mat_t A, long prec)
Computes the determinant of the matrix, using Gaussian elimination
with partial pivoting. If at some point an invertible pivot element
cannot be found, the elimination is stopped and the magnitude of the
determinant of the remaining submatrix is bounded using
Hadamard's inequality.

344
doc/source/fmprb_poly.rst Normal file
View file

@ -0,0 +1,344 @@
fmprb_poly.h -- polynomials of real number balls
===============================================================================
An *fmprb_poly_t* represents a polynomial over the real numbers,
implemented as an array of coefficients of type *fmprb_struct*.
Most functions are provided in two versions: an underscore method which
operates directly on pre-allocated arrays of coefficients and generally
has some restrictions (such as requiring the lengths to be nonzero
and not supporting aliasing of the input and output arrays),
and a non-underscore method which performs automatic memory
management and handles degenerate cases.
Types, macros and constants
-------------------------------------------------------------------------------
.. type:: fmprb_poly_struct
.. type:: fmprb_poly_t
Contains a pointer to an array of coefficients (coeffs), the used
length (length), and the allocated size of the array (alloc).
An *fmprb_poly_t* is defined as an array of length one of type
*fmprb_poly_struct*, permitting an *fmprb_poly_t* to
be passed by reference.
Memory management
-------------------------------------------------------------------------------
.. function:: void fmprb_poly_init(fmprb_poly_t poly)
Initializes the polynomial for use, setting it to the zero polynomial.
.. function:: void fmprb_poly_clear(fmprb_poly_t poly)
Clears the polynomial, deallocating all coefficients and the
coefficient array.
.. function:: void fmprb_poly_fit_length(fmprb_poly_t poly, long len)
Makes sures that the coefficient array of the polynomial contains at
least *len* initialized coefficients.
.. function:: void _fmprb_poly_set_length(fmprb_poly_t poly, long len)
Directly changes the length of the polynomial, without allocating or
deallocating coefficients. The value shold not exceed the allocation length.
.. function:: void _fmprb_poly_normalise(fmprb_poly_t poly)
Strips any trailing coefficients which are identical to zero.
.. function:: void fmprb_poly_zero(fmprb_poly_t poly)
.. function:: void fmprb_poly_one(fmprb_poly_t poly)
Sets *poly* to the constant 0 respectively 1.
Conversions
-------------------------------------------------------------------------------
.. function:: void fmprb_poly_set_fmpz_poly(fmprb_poly_t poly, const fmpz_poly_t src, long prec)
.. function:: void fmprb_poly_set_fmpq_poly(fmprb_poly_t poly, const fmpq_poly_t src, long prec)
Sets *poly* to *src*, rounding the coefficients to *prec* bits.
Input and output
-------------------------------------------------------------------------------
.. function:: void fmprb_poly_printd(const fmprb_poly_t poly, long digits)
Prints the polynomial as an array of coefficients, printing each
coefficient using *fmprb_printd*.
Comparisons
-------------------------------------------------------------------------------
.. function:: int fmprb_poly_contains_fmpq_poly(const fmprb_poly_t poly1, const fmpq_poly_t poly2)
Returns nonzero iff *poly1* contains *poly2*.
.. function:: int fmprb_poly_equal(const fmprb_t A, const fmprb_t B)
Returns nonzero iff *A* and *B* are equal as polynomial balls, i.e. all
coefficients have equal midpoint and radius.
Arithmetic
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_add(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
Sets *{C, max(lenA, lenB)}* to the sum of *{A, lenA}* and *{B, lenB}*.
Allows aliasing of the input and output operands.
.. function:: void fmprb_poly_add(fmprb_poly_t C, const fmprb_poly_t A, const fmprb_poly_t B, long prec)
Sets *C* to the sum of *A* and *B*.
.. function:: void _fmprb_poly_mullow(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long n, long prec)
Sets *{C, n}* to the product of *{A, lenA}* and *{B, lenB}*, truncated to
length *n*. The output is not allowed to be aliased with either of the
inputs. We require `\mathrm{lenA} \ge \mathrm{lenB} > 0`,
`n > 0`, `\mathrm{lenA} + \mathrm{lenB} - 1 \ge n`.
As currently implemented, this function puts each input polynomial on
a common exponent, truncates to prec bits, and multiplies exactly over
the integers. The output error is computed by cross-multiplying the
max norms.
.. function:: void fmprb_poly_mullow(fmprb_poly_t C, const fmprb_poly_t A,
const fmprb_poly_t B, long n, long prec)
Sets *C* to the product of *A* and *B*, truncated to length *n*.
.. function:: void _fmprb_poly_mul(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
Sets *{C, n}* to the product of *{A, lenA}* and *{B, lenB}*, truncated to
length *n*. The output is not allowed to be aliased with either of the
inputs. We require $\mathrm{lenA} \ge \mathrm{lenB} > 0$, $n > 0$.
This function currently calls *_fmprb_poly_mullow*.
.. function:: void fmprb_poly_mul(fmprb_poly_t C, const fmprb_poly_t A, const fmprb_poly_t B, long prec)
Sets *C* to the product of *A* and *B*.
.. function:: void _fmprb_poly_inv_series(fmprb_struct * Qinv, const fmprb_struct * Q, long len, long prec)
Sets *{Qinv, len}* to the power series inverse of *{Q, len}*. Uses Newton iteration.
.. function:: void fmprb_poly_inv_series(fmprb_poly_t Qinv, const fmprb_poly_t Q, long n, long prec)
Sets *Qinv* to the power series inverse of *Q*.
.. function:: void _fmprb_poly_div(fmprb_struct * Q, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
.. function:: void _fmprb_poly_rem(fmprb_struct * R, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
.. function:: void _fmprb_poly_divrem(fmprb_struct * Q, fmprb_struct * R, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
.. function:: void fmprb_poly_divrem(fmprb_poly_t Q, fmprb_poly_t R, const fmprb_poly_t A, const fmprb_poly_t B, long prec)
Performs polynomial division with remainder, computing a quotient `Q` and
a remainder `R` such that `A = BQ + R`. The leading coefficient of `B` must
not contain zero. The implementation reverses the inputs and performs
power series division.
.. function:: void _fmprb_poly_div_root(fmprb_struct * Q, fmprb_t R, const fmprb_struct * A, long len, const fmprb_t c, long prec)
Divides `A` by the polynomial `x - c`, computing the quotient `Q` as well
as the remainder `R = f(c)`.
Product trees
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_product_roots(fmprb_struct * poly, const fmprb_struct * xs, long n, long prec)
.. function:: void fmprb_poly_product_roots(fmprb_poly_t poly, fmprb_struct * xs, long n, long prec)
Generates the polynomial `(x-x_0)(x-x_1)\cdots(x-x_{n-1})`.
.. function:: fmprb_struct ** _fmprb_poly_tree_alloc(long len)
Returns an initialized data structured capable of representing a
remainder tree (product tree) of *len* roots.
.. function:: void _fmprb_poly_tree_free(fmprb_struct ** tree, long len)
Deallocates a tree structure as allocated using *_fmprb_poly_tree_alloc*.
.. function:: void _fmprb_poly_tree_build(fmprb_struct ** tree, const fmprb_struct * roots, long len, long prec)
Constructs a product tree from a given array of *len* roots. The tree
structure must be pre-allocated to the specified length using
*_fmprb_poly_tree_alloc*.
Composition
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_compose_horner(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long prec)
.. function:: void fmprb_poly_compose_horner(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
.. function:: void _fmprb_poly_compose_divconquer(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long prec);
.. function:: void fmprb_poly_compose_divconquer(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
.. function:: void _fmprb_poly_compose(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long prec)
.. function:: void fmprb_poly_compose(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long prec)
Sets *res* to the composition `h(x) = f(g(x))` where `f` is given by
*poly1* and `g` is given by *poly2*, respectively using Horner's rule,
divide-and-conquer, and an automatic choice between the two algorithms.
The underscore methods do not support aliasing of the output
with either input polynomial.
.. function:: void _fmprb_poly_compose_series_horner(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long n, long prec)
.. function:: void fmprb_poly_compose_series_horner(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long n, long prec)
.. function:: void _fmprb_poly_compose_series_brent_kung(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long n, long prec)
.. function:: void fmprb_poly_compose_series_brent_kung(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long n, long prec)
.. function:: void _fmprb_poly_compose_series(fmprb_struct * res, const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2, long n, long prec)
.. function:: void fmprb_poly_compose_series(fmprb_poly_t res, const fmprb_poly_t poly1, const fmprb_poly_t poly2, long n, long prec)
Sets *res* to the power series composition `h(x) = f(g(x))` truncated
to order `O(x^n)` where `f` is given by *poly1* and `g` is given by *poly2*,
respectively using Horner's rule, the Brent-Kung baby step-giant step
algorithm, and an automatic choice between the two algorithms.
We require that the constant term in `g(x)` is exactly zero.
The underscore methods do not support aliasing of the output
with either input polynomial.
Evaluation and interpolation
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_evaluate(fmprb_t res, const fmprb_struct * f, long len, const fmprb_t a, long prec)
.. function:: void fmprb_poly_evaluate(fmprb_t res, const fmprb_poly_t f, const fmprb_t a, long prec)
Sets res to `f(a)`, evaluated using Horner's rule.
.. function:: void _fmprb_poly_evaluate_vec_iter(fmprb_struct * ys, const fmprb_struct * poly, long plen, const fmprb_struct * xs, long n, long prec)
.. function:: void fmprb_poly_evaluate_vec_iter(fmprb_struct * ys, const fmprb_poly_t poly, const fmprb_struct * xs, long n, long prec)
Evaluates the polynomial simultaneously at *n* given points, calling
*_fmprb_poly_evaluate* repeatedly.
.. function:: void _fmprb_poly_evaluate_vec_fast_precomp(fmprb_struct * vs, const fmprb_struct * poly, long plen, fmprb_struct ** tree, long len, long prec)
.. function:: void _fmprb_poly_evaluate_vec_fast(fmprb_struct * ys, const fmprb_struct * poly, long plen, const fmprb_struct * xs, long n, long prec)
.. function:: void fmprb_poly_evaluate_vec_fast(fmprb_struct * ys, const fmprb_poly_t poly, const fmprb_struct * xs, long n, long prec)
Evaluates the polynomial simultaneously at *n* given points, using
fast multipoint evaluation.
.. function:: void _fmprb_poly_interpolate_newton(fmprb_struct * poly, const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
.. function:: void fmprb_poly_interpolate_newton(fmprb_poly_t poly, const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most *n* that interpolates
the given *x* and *y* values. This implementation first interpolates in the
Newton basis and then converts back to the monomial basis.
.. function:: void _fmprb_poly_interpolate_barycentric(fmprb_struct * poly, const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
.. function:: void fmprb_poly_interpolate_barycentric(fmprb_poly_t poly, const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most *n* that interpolates
the given *x* and *y* values. This implementation uses the barycentric
form of Lagrange interpolation.
.. function:: void _fmprb_poly_interpolation_weights(fmprb_struct * w, fmprb_struct ** tree, long len, long prec)
.. function:: void _fmprb_poly_interpolate_fast_precomp(fmprb_struct * poly, const fmprb_struct * ys, fmprb_struct ** tree, const fmprb_struct * weights, long len, long prec)
.. function:: void _fmprb_poly_interpolate_fast(fmprb_struct * poly, const fmprb_struct * xs, const fmprb_struct * ys, long len, long prec)
.. function:: void fmprb_poly_interpolate_fast(fmprb_poly_t poly, const fmprb_struct * xs, const fmprb_struct * ys, long n, long prec)
Recovers the unique polynomial of length at most *n* that interpolates
the given *x* and *y* values, using fast Lagrange interpolation.
The precomp function takes a precomputed product tree over the
*x* values and a vector of interpolation weights as additional inputs.
Differentiation
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_derivative(fmprb_struct * res, const fmprb_struct * poly, long len, long prec)
Sets *{res, len - 1}* to the derivative of *{poly, len}*.
Allows aliasing of the input and output.
.. function:: void fmprb_poly_derivative(fmprb_poly_t res, const fmprb_poly_t poly, long prec)
Sets *res* to the derivative of *poly*.
.. function:: void _fmprb_poly_integral(fmprb_struct * res, const fmprb_struct * poly, long len, long prec)
Sets *{res, len}* to the integral of *{poly, len - 1}*.
Allows aliasing of the input and output.
.. function:: void fmprb_poly_integral(fmprb_poly_t res, const fmprb_poly_t poly, long prec)
Sets *res* to the integral of *poly*.
Special functions
-------------------------------------------------------------------------------
.. function:: void _fmprb_poly_log_series(fmprb_struct * f, fmprb_struct * h, long n, long prec)
.. function:: void fmprb_poly_log_series(fmprb_poly_t f, const fmprb_poly_t h, long n, long prec)
Sets `f` to the power series logarithm of `h`, truncated to length `n`.
Uses the formula `\log f = \int f' / f`, adding the logarithm of the
constant term in `h` as the constant of integration.
The underscore method does not support aliasing of the input and output
arrays.
.. function:: void _fmprb_poly_exp_series_basecase(fmprb_struct * f, const fmprb_struct * h, long hlen, long n, long prec)
.. function:: void fmprb_poly_exp_series_basecase(fmprb_poly_t f, const fmprb_poly_t h, long n, long prec)
.. function:: void _fmprb_poly_exp_series(fmprb_struct * f, const fmprb_struct * h, long hlen, long n, long prec)
.. function:: void fmprb_poly_exp_series(fmprb_poly_t f, const fmprb_poly_t h, long n, long prec)
Sets `f` to the power series exponential of `h`, truncated to length `n`.
The basecase version uses a simple recurrence for the coefficients,
requiring `O(nm)` operations where `m` is the length of `h`.
The main implementation uses Newton iteration, starting from a small
number of terms given by the basecase algorithm. The complexity
is `O(M(n))`. Redundant operations in the Newton iteration are
avoided by using the scheme described in [HZ2004]_.
The underscore methods support aliasing and allow the input to be
shorter than the output, but require the lengths to be nonzero.
.. function:: void fmprb_poly_log_gamma_series(fmprb_poly_t f, long n, long prec)
Sets `f` to the series expansion of `\log(\Gamma(1-x))`, truncated to
length `n`.

26
doc/source/history.rst Normal file
View file

@ -0,0 +1,26 @@
History and changes
===============================================================================
For more details, view the detailed commit log
in the git repository https://github.com/fredrik-johansson/arb
* 2012-09-29 - version 0.2
* code for computing the gamma function (Karatsuba, Stirling's series)
* rising factorials
* fast exp_series using Newton iteration
* improved multiplication of small polynomials by using classical multiplication
* implemented error propagation for square roots
* polynomial division (Newton-based)
* polynomial evaluation (Horner) and composition (divide-and-conquer)
* product trees, fast multipoint evaluation and interpolation (various algorithms)
* power series composition (Horner, Brent-Kung)
* added the fmprb_mat module for matrices of balls of real numbers
* matrix multiplication
* interval-aware LU decomposition, solving, inverse and determinant
* many helper functions and small bugfixes
* 2012-09-14 - version 0.1
* 2012-08-05 - began simplified rewrite
* 2012-04-05 - experimental ball and polynomial code

40
doc/source/index.rst Normal file
View file

@ -0,0 +1,40 @@
Welcome to Arb's documentation!
===============================
Arb is a C library for arbitrary-precision
floating-point ball arithmetic, developed by
Fredrik Johansson (fredrik.johansson@gmail.com).
The git repository is https://github.com/fredrik-johansson/arb/
A `PDF version <http://fredrikj.net/arb/arb.pdf>`_ of this documentation
is available.
General information
::::::::::::::::::::
.. toctree::
:maxdepth: 2
overview.rst
setup.rst
history.rst
Module documentation
::::::::::::::::::::
.. toctree::
:maxdepth: 2
fmpr.rst
fmprb.rst
fmprb_poly.rst
fmprb_mat.rst
Credits and references
::::::::::::::::::::::::
.. toctree::
:maxdepth: 2
credits.rst

50
doc/source/overview.rst Normal file
View file

@ -0,0 +1,50 @@
Feature overview
===============================================================================
Ball arithmetic, also known as mid-rad interval arithmetic, is an
extension of floating-point arithmetic in which an error bound is
attached to each variable. This allows doing rigorous computations
over the real numbers, while avoiding the overhead of
traditional (inf-sup) interval arithmetic at high precision,
and eliminating much of the need for time-consuming
and bug-prone manual error analysis associated with
standard floating-point arithmetic. (See for example [Hoe2009]_.)
At the moment, Arb contains:
* A module (fmpr) for correctly rounded arbitrary-precision
floating-point arithmetic. Arb numbers have a few special features, such
as arbitrary-size exponents (useful for combinatorics and asymptotics) and
dynamic allocation (facilitating implementation of hybrid
integer/floating-point and mixed-precision algorithms).
* A module (fmprb) for real ball arithmetic, where a ball is
implemented as a pair of fmpr numbers.
* Functions for fast high-precision computation of some mathematical constants,
based on ball arithmetic.
* A module (fmprb_poly) for polynomials or power series over the real numbers,
implemented using balls as coefficients, with fast polynomial multiplication.
* A rudimentary module (fmprb_mat) for matrices over the real numbers,
implemented using balls as coefficients.
Planned features include: transcendental functions and more extensive
polynomial and matrix functionality, as well as support for complex numbers.
Arb uses `MPIR <http://mpir.org>`_ and `FLINT <http://flintlib.org/>`_
for the underlying integer arithmetic and other functions.
The code conventions borrow from FLINT, and the project might get
merged back into FLINT when the code stabilizes in the future.
Arb also uses `MPFR <http://mpfr.org/>`_, mainly for testing purposes
and fallback code.
The current version of Arb implements most of its floating-point arithmetic
naively using high-level FLINT types. The speed at low precision is far from
optimal, and the memory management can sometimes be wasteful. The internals
will be rewritten in the future to fix the inefficiencies,
which eventually should make Arb ball arithmetic about as fast as mpz or mpfr arithmetic at any precision.
**Warning**: as this is an early version, any part of the interface is
subject to change! Also be aware that there are known and unknown bugs.

44
doc/source/setup.rst Normal file
View file

@ -0,0 +1,44 @@
Installation and usage basics
===============================================================================
Arb has the following dependencies:
* MPIR (http://www.mpir.org)
* MPFR (http://www.mpfr.org)
* FLINT (http://www.flintlib.org)
Currently a source checkout of FLINT from
https://github.com/fredrik-johansson/flint2 is required.
To compile, test and install Arb from source, do::
./configure <options>
make
make check
make install
If MPIR, MPFR or FLINT is installed in some other location than
the default path /usr/local, pass the
flag --with-mpir=... --with-mpfr=... or --with-flint=... with
the correct path to configure (type ./configure --help to show
more options).
Here is a simple sample program to get started using Arb:
.. code-block:: c
#include "fmprb.h"
int main()
{
fmprb_t x;
fmprb_init(x);
fmprb_const_pi(x, 50 * 3.33);
fmprb_printd(x, 50); printf("\n");
fmprb_clear(x);
}
The output should be something like the following::
3.1415926535897932384626433832795028841971693993751 +/- 4.2764e-50