museoa/ries: RIES - Find Algebraic Equations, Given Their Solution. Retrieved from https://mrob.com/pub/ries/src/ries.c.txt. - Codeberg.org

 ries.c
 RIES -- Find Algebraic Equations, Given Their Solution
 Copyright (C) 2000-2026 Robert P. Munafo
 This is the 2026 May 06 version of "ries.c"
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 If you got ries.c from the website mrob.com, the GNU General Public
 License may be retrieved from mrob.com/ries/COPYING.txt
 You may also find a copy of the license at www.gnu.org/licenses/
 The remainder of this large header comment is broken up into sections
 titled: HOW TO BUILD, UNFINISHED WORK, DETAILED NOTES ON ALGORITHM,
 REVISION HISTORY.
HOW TO BUILD:
 1. Boot your favorite UNIX-compatible computer (Step 0: install Linux,
 because Linux rules! :-)
 2. Make a new directory, put this file there.
 3. Compile it with the following command:
 gcc ries.c -lm -o ries
 If the compilation fails and reports errors like "undefined reference",
 try moving the pieces around: "gcc -o ries ries.c -lm". Using the order
 "gcc -o ries -lm ries.c" is known to NOT work with recent versions
 of GCC.
 Try other flags for optimization if you wish (like: -m64, -O2 or -O3)
 Note that -ffast-math does *not* work; it prevents IEEE 754 compliance
 and breaks a few of the RIES algorithms.)
 4. Run ries and give it a number, for example:
 ries 2.5063
BUILD OPTIONS:
In the compile line you may add one or more of these options:
 -DRIES_WANT_LDBL
 Use this option to 'ask' RIES to use the 'long double' floating-point
 data type for all calculations. It will try to determine (via other
 predefined flags provided by the compiler) whether it is available
 and if so, will use it. This gives a few extra digits of precision on
 most Intel-based systems, and gives about 30 digits on PowerPC systems.
 On some systems (notably Cygwin) you'll get errors.
 -DRIES_USE_SA_M64
 Use this option to make RIES use standalone transcendental functions.
 These functions are in the separate source file "msal_math64.c", which
 should be available in the same place you found this source file.
 -DRIES_MAX_EXPRESSION_LENGTH=29
 Use a maximum expression length of 29 symbols, rather than the default.
 Longer expressions might be needed if you're using the --one-sided
 option a lot, but they also increase the amount of memory RIES uses
 while doing long searches.
BUILDING IN MICROSOFT VISUAL C++
If you want to build RIES in MS Visual Studio, start at the RIES
website (mrob.com/ries) and follow the "Source code" link. Download
and save the source file "ries-for-windows.c". Read its header comment
for further instructions.
RIES INSTRUCTIONS (MANUAL)
Instructions for actually using RIES are in the manual, which
should be in the same place you found this file, if not look on the
web at mrob.com/ries
The manual source is in nroff format (ries.1), and should also be
available in Postscript, PDF, and plain ASCII text. To use the nroff
version, copy it to the proper place (probably in /usr/share) e.g. :
 cp ries.1 /usr/share/man/man1
(substitute appropriate local manpage directory for your OS) then type
"man ries".
UNFINISHED WORK (TTD)
(Listed more or less in the order that I want to look at them, and
recent ones have date tags. Most of the undated notes are from prior
to December 2011)
20220318
 Add a msal_math128.c that uses the GCC __float128 type (test for
__FLT128_MANT_DIG__). Figure out how to reconcile it with library
routines (if any) e.g. sqrt can call library sqrt for the first
approximation. Perhaps also make a msal_math80.c to cover the case
when ldbl is 80-bit Intel native.
20160423
 Add --surprise-me option, which generates a target value according
to a suitably-distributed random number generator e.g.
e^(K*erf(rand(0..1))) where an approximation of the error function
would be acceptable (see en.wikipedia.org/wiki/Error_function, section
"Approximation with elementary functions")
 } else if (strcmp(pa_this_arg, "--surprise-me") == 0) {
 ries_val t;
 g_surprise_me = B_TRUE;
20150118
 There is an error in the manual: "... To exit on a match within
some ``epsilon'', use --max-match-distance with a nonzero epsilon; to
reject inexact matches use --max-match-distance ..." note that the
first behaviour does not happen (but would be a nice feature to have).
The first behaviour is accomplished by the option "-n1". Check
the rest of the manual for similar errors and make sure "-n1" and
--{min|max}-match-distance cross-reference each other.
2013.0314 Look into adding two more classes of numbers: -e for
"elementary" and something like -t for "transcendental". The
uncertainty has to do with which definitions I wish to use.
 -e for "elementary" fills the rather wide conceptual gap between -a
(algebraic) and -l (Liouvillian). I would like to have a class of
numbers where exponentiation is unrestricted, so sqrt(2)^sqrt(2) would
be allowed. That is kind of like what we'd get now with "-a
--any-exponents", but I don't want to allow x in exponents because I
consider the root of "x^x=2" to be non-elementary. So we need a
"--no-x-in-exponents" option to do -e the way I envision it. Since
trig functions are linked to exponents (Euler's formula; complex
exponential function) it seems to make sense to have a
"--no-x-in-trig" setting as well (which conveniently also makes
try.solve's job easier). I should allow e and the exponential function
but with no x in the exponent; logarithms should also be allowed, but
with no x inside a logarithm (otherwise we get things like "x*ln(x)=3"
which is just another form of "x*e^x=e^3"). This only needs one
setting, that is, --no-x-in-exponents also implies no x in arguments
to [E], [l] or in either argument of [L].
 -t would allow everything, and would also enable the W function if
it is available. The only difference versus -EW is that -t would not
generate an error when W is not available, but would just silently
proceed to find W-less solutions.
2013.0314: RIES has trouble finding an equation for
1.3566631192732151980, which is the root of 2^x+3^x=7. The best I
could do was "ries -p 1.3566631192732151980 -Sx+-23^5 -l5" which finds
the awkward "x-(x+2^x+3^x) = (2+2)^(2-3^3)-(2+5)". RIES suffers from
the LHS always having to start with [x). I could address this by
removing the LHS-RHS distinction when calling gen_forms, and find a
different way to regulate the balance between LHS's and RHS's in the
database. A change like this might be easier (or harder) is it is done
after (or before) I generalize the handling of restricted symbolsets.
It also clearly affects the "x on both sides of the equation" change.
2013.0202 --symbol.names is nice, but I need to add definable "begin"
and "end" strings for various things e.g. superscripts, and start
looking at user-redefinable fixity and precedence for things like
mapping [abv] -> "Surd[a,b]" for Mathematica.
2012.0103: One or more custom defined constants.
 To support this and a future defined-functions feature, each
user-defined thing needs to have a symbol auto-generated for it out of
the set of as-yet-unused byte values. Preferably this auto-generated
symbol will be a unique and yet-unused ASCII letter, but it might have
to use more than one letter. That means that (in addition to mapping
functions to turn names into symbols, and additional fields in the
sym_attr_block structure), I also need a new low-level output
formatter for -F0 and some of the debug_xxx output to replace the
simple printf("%s", expr) that I use now..
 This formatter should turn non-standard symbols into '(Nm)' where
"Nm" is an abbreviated version of the user's supplied name. This
symbol can also be used in the -O, -N, and -S options with parentheses
(which have no purpose in that context), so e.g. '-Op(Eg)' would allow
pi and the user-defined symbol 'Eg' to be used once each. To avoid
ambiguity the symbol might have to be auto-generated, and to avoid
driving users nuts I can give them hints on how to avoid that.
 Constants and (eventually) functions can share a common FORTH-like
syntax. Note the optional weight and abbreviation fields just before
the full name:
 # x e^-(x^2) is the inverse of Gosper's "Dilbert Lambda Function"
 --define : InverseDilbertLambda ( x -- x e^-(x^2) )
 dup dup* neg exp *
 ;
 --define : XeX ( x -- x*e^x ) dup exp * ;
 --define : Eg:EulerGamma # seft-a (constant)
 ( -- The Euler-Mascheroni constant, 0.57721... )
 ( --value-type TRAN )
 # 50 digits for when RIES goes to higher precision
 0.57721566490153286060651209008240243104215933593992
 ;
 --define : 14:H:Hypot # seft-c (two-argument function), weight 14
 ( a b -- sqrt(a^2+b^2) )
 dup* # ( a b^2 )
 swap # ( b^2 a )
 dup* + sqrt
 ;
 Note that the FORTH comment delimiter is parentheses, which fits
nicely with the fact that parentheses are unneeded in postfix
expressions.
 In the first implementation of functions, you cannot include literal
constants in the function, but instead have to --define the constant
first with its own name.
 2014.1122: Everything in comments would be ignored except a token
starting with '--', like the ( --value-type TRAN ) example above. This
allows me to extend the FORTH syntax to deliver optional metadata. If
a constant has no explicitly given --value-type, it would be guessed
using guess_valtype() as is currently done with the target. Advanced
users who are making collections of constants in -p files would define
a --value-type for all.
 2012.0613: I might also want to allow immediate constants in
--eval-expression strings (and thus, in eval()) delimited by parentheses
or whitespace. If using parentheses, "ries --eval-expression '(27)q'"
would produce similar output to "ries 27 --eval-expression xq" or
"ries --eval-expression '33^q'". This improvement would be done along
with a loosening of the just-mentioned restriction on using immediates
within functions.
2013.0228 Add DUP, SWAP and OVER operators. These require a bit of
modification to the stack depth and undo handling, but they are
similar enough to seft-a and seft-b that it shouldn't be too hard.
With a default complexity about the same as a small integer, they
would allow duplication of common subexpressions, which seems to
happen fairly often in real formulas.
2013.0314 More creative use of attribute tags (TAG_INT, TAG_RAT, etc):
 --only-integer-exponents: exponents must be integer if the
argument does not contain x, and roots must be integer if the
argument *does* contain x, and all others are disallowed.
 2013.0309:
 --unit-fractions option: reciprocal only if argument is integer
(possibly useful for Egyptian fractions work)
2012.1222 Allow all odd integer roots of a negative argument (3
currently allowed, but better if any odd integer were detected, tags
should help with this; requires extra wizardry in try.solve); add an
option to display "3,/(...)" as "cbrt(...)" (Unicode handled separately
with a .ries profile)
2016.0131 The "unicode.ries" profile is disappointing for a few reasons:
 * The legend at the bottom of RIES' output still uses the standard
 ASCII symbols
 * sqrt(2) becomes "/(2)" (where / is the radical sign) but the parens
 shouldn't be there if the argument is only one character long
 * Nth root, e.g. 3,/7 for cube root, should use a "3" superscript or
 the Unicode cobe-root symbol (U+221B), and other special cases
For all of these we need some sort of user-definable cascading list of
transformation rules.
2013.0205 Output format option that prints any integer-valued
subexpression as an integer, so that "ries -s -l3 351.36306009596398"
can give "8 sqrt(1929)" instead of "8 sqrt((5-7^2)^2-7)". Attribute
tag should help with this, but infix conversion will need to call
eval() on subexpressions to figure out when a given subexpression can
be substituted with an integer (and recursive calls to infix
conversion pruned).
2012.1222: In --one-sided mode we could take multiple targets, enter
them all in the database as [x], and use the existing algorithm which
would thereby compare every match to each target. This requires a
significant change to exec() in that it must get its exec_x from a
(new) field in the pe, since the "x" value differs from one expression
to another.
 For --one-sided mode (with a single target or with multiple targets)
we really don't need a database at all: every generated expression can
just be compared to x on-the-fly. This would save memory and run all
in cache, possibly much faster, and allow for a simple multi-threaded
implementation. This can be combined with the multiple targets
improvement (in which case there would be a small database containing
the [x] expression for each target).
2013.0305: A similar idea to the "multiple targets" idea is a
"correlation search" operating mode. In this mode there are just two
targets T and U, and all expressions contain one or the other; any
match must have a T expression on the LHS and a U expression on the
RHS. This is like one-sided mode except that everything is an LHS, and
it's like the planned "x on both sides of the equation" mode in that
matches need to deal with both sides of the equation having a
derivative. One big difference is that the derivative with respect to
T might be independent of the derivative with respect to U (or they
might be partly correlated, or anti-correlated). Thus, each reported
match would report the "delta" as a vector (with direction and
magnitude) encoding how far and in what direction each of the two
targets T and U need to be altered to make that particular equation
match.
2013.0307: If exiting because of a parameter error, display a
traceback with the name(s) and character offsets of any include files
(Line numbers would be nice too, but might be difficult).
2012.1207: Review implementation of --try-solve-for-x:
 * Inverse trig functions are multi-branched. (If they give 7.012913
we might find "sin(x) = 2/3", when solving that for x we must recognize
it and solve to "x = arcsin(2/3)+2 pi" rather than just "x = arcsin(2/3)").
 * Add new options --LHS-addsym, --LHS-onlysyms, --LHS-onesym, and
--LHS-nosyms specify the symbol-set for LHS generation (and likewise
for RHS). This is a refactoring, to prepare for the next changes.
 * The -S and -N options need to have a different meaning when solving
for x. sym.attrs[i].sa_alwd (FKA sym_allowed) needs to be split into
LHS and RHS parts (to provide for the fact that, when solving, some
operators become their inverses: thus the option -SE should set
LHS.allowed['l'] and RHS.allowed['E'] to nonzero values).
 * To provide compatibility, the -S, -E, -O, and -N options work the
old way unless followed by an '=' character. For example, '-N=q' will
forbid sqrt in RHS and forbid x^2 in LHS.
 * Likewise, the symbol-related variables, like a_minw etc. need to be
duplicated so we have one set for LHS and another for RHS.
 * The sym_allowed handling needs to be even more subtle for LHS
expressions: In the expression [xl3^23/-], the [23/] will move
unchanged to the RHS, whereas the [l], [3^] and [-] will get changed
into [E], [3v] and [+] respectively on the RHS. The way to handle this
is for ge.2 to look at whether the argument(s) contain x (by testing
dx!=0) and can also take the current stack pointer into account when
deciding whether to add a symbol. For seft-b symbols like [l], if the
value on the stack contains x (or if the SP is precisely 1) then the
LHS sym_allowed array should be used, but if the SP is bigger than 1
the RHS sym_allowed should be used. Likewise the seft-c symbols should
check if the SP is precisely 2. This also implies that in
--solve-for-x mode the gf_1 and ge.1 complexity limits need to be
based on the union of LHS and RHS symbolsets.
 * Dealing with extra x's in the equation. Options include:
 - Use of -s implies or requires -Ox
 - Use of -s implies or requires the "X only on LHS" mode (opposite of
allowing x on both sides of an equation)
 - Continue with the current practice of simply pushing any extra
x's over to the RHS.
 * Dealing with trig functions. The options are:
 - Adding inverse trig functions, which would clutter the search
space even more.
 - Use of -s implicitly disables trig functions only on LHS (you can
explicitly enable them with an --LHS-syms option)
 - Use of -s implicitly disables trig functions, and/or implicitly
enables their inverses (which would be disabled by default)
 - Use of -s shifts the weights of trig functions to make them less
dense in the search space
 Perhaps these different options could be selected via an extra letter
after the -s
 * The -O option, or any similar option setting a specific non-zero
limit on the number of a particular symbol, cannot be implemented
efficiently. (Consider what happens if all expressions are inserted
into the same tree: then for any given expression, a large fraction of
the other expressions in the tree cannot be used to make a valid
solution: If we get "ln(x)=e^2" and solve for x, and they gave the
option -OE, that solution is not valid). To resolve this, the best
solution is probably to make "solve for x" mutually exclusive to the
-O option. As a sole exception, -Ox can still be implemented, so long
as all reported results involve just a single LHS and an RHS.
 * 2013.0215: One-sided simplification/reduction rules similar to
those already in cv.simplify. One simple example is [ss]->[4^]. A lot
of these are already implicit in the pruning rules, so I can look in
the comments next to the add_rule() calls to get the list of
simplifications. Virtually all will be simple string substitutions.
The unsimplified forms show up only when using -s.
2012.0503 (partly done on 20120505): Prior to 20120505, "ries -l-2 -i
143" gave the bizarre answer "(x-2)+3 = (3*4)^2". The "(x-2)+3"
(instead of x+1) results from the k_sig_loss value 0.01: Because 1 is
less than 143*k_sig_loss, RIES doesn't want to perform the addition
x+1. Similarly, "ries 7775" fails to find "x+1=6^5" because it doesn't
want to add 1 (or any integer) to x.
 The narrower, safer fix (implemented on 2012.0505) is to attempt
"(x+1)-1" and see if it is precisely equal to x. If so, then we can
assert that no significance loss has actually occurred, and permit an
exception to the k_sig_loss rule for addition/subtraction. Other
similar "conservative" tests might be possible, but add/subtract is
the big one.
 A broader fix to investigate is to always allow f(x)+K whenever the
variable subexpression f(x) is larger in magnitude than the constant
subexpression K. That will require testing with #qualify#.
2012.0511 (partly done on 20120720): The "four 4's problem" and many
others like it are small enough to implement with --one.sided and/or
--numeric.anagram. RIES's handling of these problems can be improved
with fairly little modification:
 * An option to not use rules like [xx/] that are meant to prune
"redundant" subexpressions like "4/4".
 * As a separate option (perhaps by giving "--numeric.anagram" vs.
"--strict-numeric-anagram"), accept an expression only if all of the
symbols have been "used up" in the expression. This will require
changing the "exhaustion timeout" detection: I can probably use a
"Still searching" message similar to the one I added to handle
--max-match.distance.
2012.0520: Presently the report.match routine has a bunch of arrays
called {r|l}_{e|f|g}scratch, and the escratch ones are being used for
two different purposes (char[] and symbol[]). The names should be
improved, and a dedicated symbol[] scratch array added. These changes
are logically combined with -Fmax, -FHTML, etc. output formats:
2013.0314: Augment --symbol.names in whatever ways are needed for a
user to implement output formats such as: raw symbols,
calculator-keys, HTML, RHTF, TeX, eqn, Maxima, etc.
 Each symbol should have a user-definable precedence, associativity
and layout (like Haskell "fixity"). If desired, an operator like *
(seft 'c', FORTH stack effect (a b -- c) ) could be redefined to be
displayed as any of the following:
 arg_order s_0 s_1 s_2
 Lisp: (times a b) a, b '(times ' ' ' ')'
 RIES: ba* b, a '' '' '*'
 infix: a x b a, b '' ' x ' ''
 HTML: a×b a, b '' '×' ''
There also need to be flags indicating when things can be left out, as
is presently done with '*' in certain cases, and indicating when
parentheses can be left out, as is presently done with a product
inside a sum.
 It is important to note that user-defined functions leaving more
than one item on the FORTH stack cannot be expressed in infix. (If the
two outputs of a single function are immediately multipled together,
what goes on each side of the multiplication sign? It's even worse
when things are done to some of the outputs before they are combined.)
2011.1226: Consolidate the (current) multiple ways that roundoff, overflow,
loss of significance, and tautology errors are handled. This involves the
constants and variables: n_ovr, p_ovr, k_sin_clip, k_min_best.match,
k.vanished_dx, k_prune.deriv, and k_sig_loss. For example,
k_sig_loss*k_min_best.match should be equal to the magnitude of the
"machine epsilon" (2^-53 in the case of 64-bit double) determined by
init.formats().
2012.0106: More importantly, we need to allow matches between an LHS and
another LHS (producing solutions with X on both sides of the
equation). The match closeness becomes |(val_r - val_l)/(deriv_l -
deriv_r)| which is the same as the current formula except for the
addition of the deriv_r term.
 * (2011.1226) Start at the bottom (exec(), eval(), newton(), etc.) and work my
way up to the top (ge.1() etc.). Rename variables and parameters (like
the "!on_rhs" passed to exec() from ge.2) so the names agree with what
they actually do ("!on_rhs" should be something like "!no_x" or
"has_x").
 * Allow more user control over how much output contains x on both
sides, from old behaviour to "all eqns have x on both sides". The
default should be somewhere in the middle; this should be implemented
in the main outer loop where it tests "if (lhs_insert > rhs_insert)".
Also, the 'solve for x' and '-Ox' options should both force the old
behaviour (the latter because of the need to have different
sym.attrs[i].sa_alwd values for LHS vs. RHS).
2011.1226: Higher precision:
 * (partly implemented on 2013.0301) On Intel targets in GCC 4.2.1,
__LDBL_MANT_DIG__ is 64, implying that I can use "long double" to get
a bit more precision. To exploit this I will need to do runtime
testing to determine what precision I am actually getting, similar to
the code in hint.c
 * (complete by 2013.0626) Continue conversion to use of long double
and runtime adjustment of epsilons and cutoffs.
 * Add string-to-float and float-to-string conversion routines based
on f107_spfg and f107_sscan in f107_o.cpp. (I think these already
handle the exponent adequately, but I'll want smarter handling of
extra whitespace when there is no exponent). Convert all instances of
double printf and scanf to use these routines (including the debug_X
printfs).
 * Use three sets of flags:
 Desired precision: RIES_WANT_F53, RIES_WANT_F64 and RIES_WANT_F107;
 also RIES_WANT_LDBL etc. for users who don't know what their 'long
 double' precision but still want to use it.
 Available types and their precisions: RIES_HAVE_LDBL_F64,
 RIES_HAVE_LDBL_F107 and RIES_HAVE_F128
 Which type to use: RIES_VAL_DBL, RIES_VAL_LDBL and RIES_VAL_F128
 * Don't bother with C++; just use __float128 if it's available and
too bad if it's not.
 * No need to use macros like MUL(src1,src2,dst) because __float128 and
'long double' have fully implemented builtins.
 * Initialize the special constants (k.vanished_dx, k_prune.deriv,
etc.) differently depending on precision option. It should be possible
to add these gradually without breaking normal precision since the
rest of the code will just ignore the lower half of any quad
variables.
Consider options for utilizing multiple processors (threads / parallel
implementation). {As of 2012.0423} I have identified three possible
paths of development, below titled "Best", "Good" and "BAD":
2012.0423: Best: Perform additional searching *without* expanding the database.
This would be done when the memory limit is reached (perhaps in
response to a user option):
 - Walk the tree, applying all possible monadic transforms to each
node, to see if the result then produces a new match to another
existing node. For these purposes a "monadic transform" consists of
appending *either* a single seft-b symbol, *or* a seft-a followed by a
seft-c. For example, if this scan found [43L] in the tree, it would
append 's' to form the expression [43Ls], compute its value (which is
1.592289...) then search the tree for this value using a bt_find()
function. It might then find [xp/], which would be a near-match in the
case of x=5. This type of search can be efficiently parallelized by
having each of N threads traverse 1/N part of the tree. (Which in turn
suggests that we should maintain population-counts at each node and
add a bt_index() routine.)
 - Greater-complexity expressions can be synthesized by appending a
short, low-complexity complete subexpression and a seft-c symbol.
For example, [43L] : [3q] : [+] => [43L3q+], which would then match
[x2-].
2012.0109: Good, but hard to implement, and very memory-intensive:
 Within each petit-cycle, have N threads
running at any one time, each constraining itself to a part of the
expression search-space, distinguished by the first few symbols in the
expression. For example, after the complexity depth has gotten high
enough, we will have identified all possible combinations for the
first 3 symbols in the expression's postfix representation. The
searchspace can then be partitioned by having ge.2 perform a CRC
hashfunction on the first 3 symbols, and determine whether it should
prune or recurse depending on the low log_2(N) bits of the hash value.
 This requires having each task build up its own binary tree, which
in turn requires a parallel binary tree merge. For example, using 4
threads, each with its own binary tree divided into quartiles:
 * Add population-counts to each tree node and add a bt_index()
routine;
 * When it is time to merge, each of the 4 trees becomes read-only;
 * Each of 4 threads then builds up a new tree consisting of the I'th
quartiles of each of the old trees;
 * Combine these 4 new trees into a single big tree by creating 2 dummy
parent nodes and 1 dummy grandfather node;
 * Deallocate the old trees.
 The main problem with this approach is that it uses twice as much
memory during the merge process, and 4 times as much memory bandwidth;
and (more crucially) there is no clear way to determine in advance
whether there will be enough free memory to perform the merge.
2014.1122: Good, and a little less hard to implement:
 Run in normal single-threaded mode until the tree is big enough to
survey the location of quartiles as in the 2012.0109 proposal. Then,
permanently split the tree into N pieces (probably by a simple
traverse-and-copy). On subsequent patit-cycles, run N threads where
each thread imposes its own values of g_min_equ_val and g_max_equ_val.
Threads will find duplicate solutions, but are not fully redundant; as
a rough guess I imagine it will have a factor of sqrt(N) of
ineffeciency, so a 4-thread system will use 2 times as much memory
to find the same answers in half as much time.
2012.0423: BAD: The following does *not* work because of the elaborate
interleaved nature of the bt.insert algorithm (see note in oldries.c
mentioning "memory performance degradation"), however it could if we
return to the older 2-tree implementation.
 - Have 2 scans running at any one time, an RHS scan and a LHS scan.
Whenever a scan completes, initiate another. Every time a scan
completes, another thread runs through the outputs of the latest RHS
and LHS.
%%% options to add:
2009.0603: Higher complexity scores for trig functions unless argument
expression contains pi or x. {This was mostly addressed by the trig
argument scale change on 2011.1229}
2012.0102: Warren Smith suggests, "I dont mind sin(1), but sin(1) IN
AN EXPONENT like 3^sin(1) is just ridiculous.", suggesting an idea
like tables of rules that match the end of a forth subexpression (1S^
in this example) and if matched, either prune it outright or add to
the complexity score. (If adding to the complexity score, the
possibility of this has to be considered by the bounds-setting and
short-circuit recursion code). After further discussion he suggested
that no subexpression should be completely excluded, but an "entropy
function" should be used to weight entire expressions based on the
likelihood they would be found in actual maths. Things like 1S^ would
get a really high weight, and the expression database would be
"exponentially" more efficient at covering the types of expressions he
is interested in.
 This is distinct from --rational-trig-args, which only allows
all-or-nothing pruning. However, Warren's idea could be used by adding
an Identity-like operator that has a symbol weight and has the effect
of adding a "blessed" tag to the value. Once blessed, the value would
then be eligible for use in an exponent.
2012.0103: --log-base option to do for ln what --trig-argument-scale
does for sine and cosine.
Ideas from 2007.0703 or earlier (all of these are in MBP's first
backup, old PMG5 archives seem to have nothing; old iMac-G4 archives
might have more info):
variations on -l that allow specifying limit of search by
precision, by time, by memory usage, or by number of equations tested.
The present -l is usually correlated with all of these but is never
equal to any of them. ('-lmem=10M', '-ltime=20m', '-ldigits=8',
'-lexprs=3e6', -leqns=1e10', etc.) --max.memory is a different because
a small -l will still cause RIES to exit before the indicated amount
of memory is used.
%%% command-line option: if no symbol '1', then 'r' shouldn't be
allowed. Likewise for '2' + '^' and 's'; '-' and 'n'; 'e' + '^' and
'E'; others? Which should be the default, the current simpler behavior or the
"more correct" behavior? The "simple" behavior is useful, as
illustrated by the 1.4142135 example in the manpage, because it
helps people find alternate ways to express the same solution.
macros (user-defined constants and functions):
 - use a separate symbolset namespace (implies macros cannot call
themselves or each other, which avoids lots of problems)
 - To use -O with a macro, include '_' in the -O symbolset
list; all symbols after '_' are macro names. -S assumes all macros
(why would you define a macro and not want to use it?); using -N
with a macro is a contradiction.
 - Macro can use any symbols, and symbols used within macro
expansions don't count against their sa.alwd (FKA sym_allowed) quotas
(this is probably easy; just need to make sure)
 - Need to implement stack-overflow detection in the metastack
routines, both for the stack and for the undo lists.
 - Test evaluator routine needs to make sure macros fit one of the
three allowed types ('a', 'b' or 'c') and that only types 'b' and 'c'
dip into the stack's current contents.
 - Type 'a' can be "precompiled" since they amount to custom
constants
 - Should I allow immediate operands (like 0.577...) to facilitate
defining constants that can't be computed any other way? How about
special operators, like 'exch'?
 - execution can probably be accomplished by having eval() call
itself recursively. It needs to handle its own error returns.
 - with many -i test cases the RHS expression totals are a lot bigger
than the LHS totals (while the insert totals are equal). There are
probably cases that come out the other way around, too. In cases like
this, the program could probably find more matches more quickly by
shifting the balance between LHS/RHS to favor the side that is
generating and inserting more expressions per unit time. (The present
implementation tests "lhs_insert > rhs_insert" which is usually the
best way to do it.) I could model this with formulas that take into
account the total amount of dead-ends, expressions, and inserts on
both sides, compared as a ratio to the equation total (and take the
1st derivative to figure out what side of the maximum you're on). The
idea is to maximize the number of full equations found per unit time
rather than keeping the number of LHS inserts equal to the number of
RHS inserts.
 - here's a way to allow printing more than one exact solution: When
inserting an LHS or RHS, if an identical item already exists in the
tree, overwrite it with the new one. This will cause at most one new,
distinct exact match per matchscan. Since this new match is of higher
aggregate complexity than the old one, it has at least a moderate
chance of being a mathematically distinct solution (-:
 - spend a lot of time looking at which expressions are generated
first, and try to improve the weights so it makes more sense. Why does
2.5063 find [x2q/]=[pq] before [xs]=[p2*]?
 - add constants: Feigenbaum, Euler's
 - add Gamma function (many notes on this below) with an option to
use "factorial" notation when displaying.
 - explore how to expand RIES to complex numbers and complex analytic
functions. Looks easy -- after all, complex data types are native in
GCC!
(Test cases are now in #qualify#)
TABLE OF FUNCTIONS
 sym stack-effect description
 0x1 -- Phantom symbol for argument-reversed version of -
 0x2 -- Phantom symbol for argument-reversed version of /
 0x4 -- Phantom symbol for argument-reversed version of ^
 ' ' -- Blank space: no operation (for --eval-expression)
 0 (-- 0) The constant 0.0 (not currently used, but will have a
 role soon as part of --numeric.anagram)
 1 (-- 1) The constant 1.0
 2 (-- 2) The constant 2.0
 3 (-- 3) The constant 3.0
 4 (-- 4) The constant 4.0
 5 (-- 5) The constant 5.0
 6 (-- 6) The constant 6.0
 7 (-- 7) The constant 7.0
 8 (-- 8) The constant 8.0
 9 (-- 9) The constant 9.0
 ! (a -- f) Factorial monad f(x) = x! = Gamma[x+1] (reserved, not
 yet implemented)
 # (n d -- x) Digit paste operator x(n,d) = 10*n+d (reserved, not yet
 implemented. This is to make --numeric.anagram more useful)
 % (a b -- m) (reserved for modulo or remainder?)
 ^ (a b -- f) Power: f(a,b) = a^b
 * (a b -- p) Multiply+ p(a,b) = a*b
 ( -- Comment delimiter; used to bracket custom symbol names;
 used as a placeholder during infix translation
 ) -- Comment delimiter; used to bracket custom symbol names;
 used as a placeholder during infix translation
 - (a b -- d) Subtract: d(a,b) = a-b
 + (a b -- s) Add: s(a,b) = a+b
 : -- Function definition start
 ; -- Function definition end
 . -- Placeholder for multiplication during infix formatting
 / (a b -- r) Divide: r(a,b) = a/b
 A (a b -- A) Arctangent A(a,b) = atan2(a,b) (reserved, not yet used)
 C (x -- c) Cosine c(x) = cos(pi x)
 E (x -- f) Exponential function f(x) = e^x
 G (x -- g) Gamma function g(x) = Gamma[x] = (x-1)! (reserved, not
 yet implemented)
 I (x -- x) Identity function x(x) = x (not used in expressions, but
 used as a placeholder during infix translation)
 L (a b -- l) Arbitrary logarithm l(a,b) = log_b(a)
 S (x -- s) Sine s(x) = sin(pi x)
 T (x -- t) Tangent t(x) = tan(pi x)
 W (x -- w) Lambert W function, e.g. w(10.0) = 1.74553...
 e (-- e) The constant 2.71828...
 f (-- f) The constant 1.61803...
 l (a -- l) Natural logarithm: l(a) = ln(a)
 n (a -- n) Negate: n(a) = -a
 p (-- p) The constant 3.14159...
 q (a -- q) Square root: q(a) = sqrt(a)
 r (a -- r) Reciprocal: r(a) = 1/a
 s (a -- s) Square: s(a) = a*a
 v (a b -- v) Root: v(a,b) = a^(1/b)
 x (-- x) The user's target number
ARCHITECTURE
 main -- Outer loop -- increment complexity and decide whether to add to
 RHS tree or LHS tree
 gen_forms -- setup first recursive level of gf_1
 gf_1 -- Add a symbol of type 'a', 'b' or 'c' and recurse if complexity
 : allows
 :.gf_1 -- Recursive levels until form is complete
 ge_1 -- setup metastack and initial level of ge.2
 ge_2 -- Generate one step of FORTH code and update metastack,
 : recurse if more symbols in the form
 :.ge_2 -- Recursive levels until form is complete
 canonval -- Try to put expression value in [1.0,2.0) (only if
 --canon-reduction option is used)
 bt_insert -- Insert calculated result into LHS or RHS tree
 check_exact_match (called for exact matches)
 report_match -- Display match without doing Newton
 check_sides -- Find closest neighbor of the opposite sidedness
 and determine if the pair sets a new record.
 check_match -- Check a single pair to see if it's a new
 solution
 cv_simplify -- Simplify equation by removing e.g. "2*"
 from both sides
 newton -- Use Newton's Method to locate ideal value of
 x for a given solution
 report_match -- Display a match that converged by
 Newton.
 try_solve -- Try to convert "expr1 = expr2" into
 "x = bigexpr"
DETAILED NOTES ON ALGORITHM
To reduce the search time, the graph theory "bidirectional search"
(also called "meet-in-the-middle") technique is used. Rather than
attempting to search for solutions of the form
 X = expression (1)
it searches for solutions of the form
 f(X) = expression (2)
where f(X) represents an expression containing X. Each side of
equation (2) is represented as a set (in memory, a ordered list,
implemented with a binary tree) of expressions and their values. Thus
there are two lists {For efficiency, these two lists are stored in a
single binary tree, so that I can check for a match after every new
insert. With a single tree, matching items end up being close to each
other in the tree}. Each list is kept in order by numerical value.
Then, the lists can easily be scanned to locate any matches. A match
consists of an element in the left-hand list whose value is close to
that of an element of the right-hand list. Keeping the lists in order
also allows duplicate expressions, like "2*pi", "pi+pi", and "pi*2",
to be ignored, because they will end up having the same value and
therefore end up being in the same spot in the list.
It is important to look for "simpler" matches before going to more
complex ones. Thus, expressions have a complexity score, computed by
assigning a certain number of points to each of the symbols in the
exression. Matches are checked in order of increasing total complexity
(which is the complexity of the left-hand side plus that of the
right-hand side)
In order to avoid the sorting problems necessary to search the
expressions in a perfectly increasing complexity order, complexity
scores are lumped into discrete finite quanta, which are actually
implemented by using small integers for each component of a complexity
score, rather than real numbers.
This does not actually cause everything to be searched in order, even
modulo the quantization errors. Some equations can be represented in
an unbalanced form with lower aggregate complexity than the
least-complex balanced form: an example is the cube root of 2: The
solutions are "x^3 = 2" and "x = 3,/2", both are unbalanced because
there is one symbol on one side of the = sign and three symbols on the
other. (",/" is the "Nth root" symbol). There are no good balanced
solutions. In such cases the unbalanced, low-complexity form will show
up later in the search than it "should".
To make generation and evaluation easy, expressions are represented in
a FORTH-like syntax with one character per symbol. Thus, "11+" is 1+1,
"2q" is sqrt(2), "ep^" is e^pi, etc. Symbols are categorized by their
stack-effect, which is abbreviated "seft" in the code.
In the actual FORTH language, the stack effect (seft) encapsulates
what a word does to the stack. It is described by a comment (some text in
parentheses) containing: the stack contents before the word is executed,
a "--" symbol, and the stack contents after the word executes.
There are three sefts in ries, 'a', 'b', and 'c'. These are the
similest and most common of a larget set of sefts found in FORTH-like
languages, of which the following are a representative sample:
 seft 0: ( -- )
 No operation (compare to b, b2 which don't change the stack level
 but do use something on the stack)
 seft a: ( -- K )
 Adds one thing (a constant) to the top of the stack
 seft a2: ( arg1 -- res1 res2 )
 Takes one value from the stack, performs an operation, and puts two
 results back on the stack. (Examples: DUP, divmod)
 seft a3: ( arg1 arg2 -- res1 res2 res3 )
 Takes two values from the stack, performs an operation, and puts three
 results back on the stack. (Example: OVER)
 seft b: ( arg -- result )
 Takes one value from the stack, performs an operation, and puts one
 result on the stack
 seft b2: ( arg1 arg2 -- res1 res2 )
 Takes two values from the stack, performs an operation, and puts two
 results back on the stack. (Examples: SWAP, polar to rectangular conversion)
 seft c: ( arg1 arg2 -- result )
 Takes two values from the stack, performs an operation, and puts one
 result on the stack (Examples: DROP, most binary operators)
The RIES symbols of each seft are:
 0: ' ' (blank space)
 a: x 1 2 p 3 e 4 5 f 6 7 8 9 (x, the digits, and the constants pi, e and phi)
 b: r s q l n S C T I (reciprocal, squared, square root, ln, negate, sine,
 cosine, tangent, identity)
 c: + - * / ^ v L (add, subtract, multiply, divide, exponent, root, logarithm)
For ease in debugging, most symbols have letters that make sense, but
not always (e.g. "f" for phi, the golden ratio; "v" for root, which
is meant to represent the v-shaped part of the standard root sign).
The only constants are the digits 1 through 9 -- no zero, and no
multidigits or decimal fractions. To get 10 you have to do "25*" or
something similar; for 1.5 you have to do "32/". A fair amount of
effort has been put into setting the symbol scores such that the score
of "25*" is only a little higher than the score for "9".
example cases for deriving the symbol weights:
 (+) ~= (*) (because both occur equally often)
 (-) > (+), but not by much
 1, 2, 3, 4, 5, ... degrade gracefully and kind of like a slide rule
 (14*) = (22*), etc. (implies (n) proportional to log(n) for 1<=n<=9)
 (33*) = (9) (implies (*) = [..] where [..] is the weight of any two symbols)
 (25*) > (9), but only by a little (no problems so far)
 [2q] ~= [5] :: therefore (2q) + [.] ~= (5)
 [x2v] = [xq] :: (2v) + [.] = (q)
 [x2^] = [xs] :: (2^) + [.] = (s)
 [99+] can be >> [55*] (because smaller numbers are more likely)
%%% complexity score for a symbol can be context sensitive, e.g. 'l'
takes a higher score the second time it is used in an equation -- so
long as its range of possible scores is within the total range for its
seft. This might be a good way to eliminate some of the nonintuitive
aspects of the current system.
Once a set of symbol scores is worked out, it can be adjusted by
adding any constant to all the values (this adds a bias for long
expressions, or a bias against long expressions. In the above
"normalized" examples there is no bias, but the shorter ones will get
generated first anyway.)
The lists start out small and grow as the program searches. On the
first pass, the left-hand list consists simply of {X}, the single-
element expression for the search value, and the right-hand list
consists only of a few common constants, for example {1, 2, *e*,
*pi*}. After a few passes the left-hand list will include things like
"xv" (sqrt(x)) and "x1+" (x+1), and the right-hand side will have
similar forms such as "2v" (sqrt(2)) and "23/" (2/3). At that point,
if X is 4/9 a match will be found between "xv" and "23/", even though
the combination of "x" and "49/" has fewer symbols overall, because
the former matching is more evenly balanced. Another way of saying
this is that, since the complexity scores of "xv" and "23/" are both
lower than the score of "49/", the "xv = 23/" matching is found before
"49/" even gets a chance to be generated.
To minimize time spent generating nonsense expressions, expressions
are generated from "forms". A "form" is a symbolic expression of
expression syntax, like "aabc" for "12q+". Each letter represents a
different type of stack operation (called "seft" above). Type "a"
pushes one item, type "b" leaves the stack with the same number of
items, and "c" leaves the stack with one less item. A form is legal if
the stack depth (the "a" count minus the "c" count) is > 0 at all
times and = 1 at the end. Thus all forms start with an "a". The number
of forms for N={1,2,3,...,8} is {1,1,2,4,9,21,51,127} (the Motzkin
numbers; Sloane's A1006). By comparison, the total number of forms
without these restrictions would be 3^(N-1):
{1,3,9,27,81,242,729,2401}. The savings is considerable. In the
left-hand list the first symbol will be an "x" in all generated
expressions. Forms are generated on the fly and never stored in
memory. As the search proceeds, longer forms are generated as needed.
For each form there is a "minimum" expression and a "maximum"
expression, as rated by complexity score. For example, "x1+" is the
minimum expression for form "aac" and "99L" (log base 9 of 9) is the
maximum. These minimum and maximum expressions are easy to find,
because each symbol has its own score and the expression score is just
the sum of the symbol score. Thus it is easy to determine, at the
beginning of each pass, which forms can generate expressions that are
within the current complexity range.
Expressions are generated from each valid form. The generation of
expressions uses a recursive backtracking algorithm, starting with the
first symbol and moving to the right. At each step, it checks to see
of the symbols we have so far still allow an expression which falls
within the complexity range -- if not, the latest symbol is dropped or
changed. It also avoids generating certain patterns of symbols which
would be of no computational value, or which are always equivalent to
a different (sometimes shorter) set of symbols. Examples of this
optimization are:
 - "aa-" for any seft a symbol would be 0, so is not generated.
 - "aa+" is the same as "a2*", so is not generated.
 - "aa*" is the same as "as", so is not generated.
 - "aa/" and "aal" are equivalent to 1, and are not generated.
 - Almost any operator after "1" is not generated ("1x", "1/", "1r",
"1^", "1v", "1l", "1L", "1v" all do nothing useful)
 - "2^" and "2v" are equivalent to "s" and "q" respectively, and are
not generated.
 - "pS" is equivalent to 0
 - "sq", "qs", "nn", "rr" all amount to nothing and are not generated.
The expressions are also evaluated while they are being generated. Any
subexpression that causes an error, such as divide-by-zero, can be
skipped along with all expressions that start with that subexpression.
For example, any expression starting with a constant followed by "nq",
such as "2nq" which means "sqrt(-2)", will be skipped.
As the search proceeds, each equation (that is, each combination of an
LHS with an RHS) is checked to see what value X would have to be for
the equation to work out perfectly, and the difference between this X
and the user's supplied number is called the "delta". An equation
becomes a "record-setter" if its delta is smaller than any delta seen
so far.
If a strict "non-fuzzy" comparison were made, due to roundoff errors
it would be possible for the same solution to be found twice. For
example, if X is near 4/9 = 0.44444.., two solutions are "x = 49/" and
"xq = 23/". One solution might be printed after the other because of
roundoff in the square-root operation. To avoid this, every time a new
record-setter is found, the criterion for another record is set to
0.999 times the new delta.
Determining the value of X for which the equation works out perfectly
is a hard problem. In theory one would have to solve the equation for
X. That involves lots of shifting and rearranging of symbols, and if
there is more than one X in the equation it can get into some rather
complicated algorithms.
Instead, this program compares the LHS and RHS values directly. (The
difference between the LHS and RHS is the "lhs/rhs diff", and varies
from the "delta"). This leads to several problems. Sometimes two sides
of the equation form a really good match only because both sides
involve something raised to a very small power (or a very big root --
same thing). An example of this type of problem is 1.017262042 =
[49sv] (the 81st root of 4) and 1.017313996 = [38sv] (the 64th root of
3). They differ by only 1 part in 20000, but only because both numbers
are very close to 1. All 81st and 64th roots of small integers are
close to 1. This is referred to as the "error margin problem" (or
"loss of significance").
There is also the problem of zero subexpressions and constant
subexpressions involving X. An example is "x1el-^ = 42sL": The
subexpression "1el-" is "1-ln(e)" which evaluates to 0, and the entire
LHS is x^(1-ln(e))" which will always be 1. The right-hand side (log
base 2^2 of 4) is also 1, so this equation works for any X. In order
for the program to work, such "solutions" have to be detected and
eliminated. This is the "tautology problem".
Finally, there is the issue of reporting what X would perfectly solve
each equation. This is a desirable feature of the program because
users will often want to find out what the "exact" value would be for
a given equation without actually doing the algebra and calculations
themselves. In cases like [xle+q] = [1e-s], where there is only one X
in the equation, this could be done by the computer. However, there are
lots of important cases like [xx^] = [52*] where an analytic solution
does not exist and the value of X has to be found by some other method,
such as a numeric method. This is the "root-finding problem"
Conveniently, all three problems (the error-margin problem, the
tautology problem, and the root-finding problem) are handled in the
same way: by calculating the *derivative* of all terms and
subexpressions on the LHS. The derivative, when multiplied by the
lhs/rhs diff, gives a very good approximation to the delta, which as
mentioned above is the amount that X would need to be changed to make
the equation a perfect match.
This solves the error margin problem because it equalizes the field --
all equations can get rated in terms of how much X has to change to
make LHS and RHS match, rather than how much the LHS would have to change.
It solves the tautology problem because, by definition, an LHS with
a zero derivative is constant with respect to X, and therefore
constitutes a tautology solution.
It also solves the root-finding problem quite conveniently. If you have
an equation like X^X = 10, and you have an approximation for X, you can
use derivatives and Newton's method to find a better approximation for
X. When RIES takes the derivative of the LHS and multiplies it by the
lhs/rhs diff, it is essentially performing one step of Newton's method.
The resulting value can be added to X to form a better approximation to
the root of the equation.
Derivatives are evaluated for LHS expressions as the expressions are
generated. In the following table, the calculation of the derivative
for each symbol is shown, based on the values of the operands (A and
B) and their derivatives (da and db). Some of these will be familiar
to anyone who has taken calculus:
 seft '0' symbols:
 ' ' no-op
 (other sefts will go here: roll, drop, dup, over, etc.)
 seft 'a' symbols:
 3 any constant 0
 x target 1
 seft 'b' symbols:
 e exponential e^A da
 l natural log da / A
 n negate - da
 s squared 2 A da
 q square root da / 2 sqrt(A)
 r reciprocal - da / A^2
 A arctan(a/b) da / (1+A^2)
 C cosine -sin(A) da
 G Gamma %%% reserved (mutex with !)
 ! factorial %%% reserved (mutex with G)
 S sine cos(A) da
 T tangent (1 + tan^2(A)) da
 W LambertW W(a)/(a(1+W(a))) da
 seft 'c' symbols:
 + plus da + db
 - minus da - db
 * times A db + B da (note "product rule" below)
 / divide B da - A db / B^2 (note "quotient rule" below)
 ^ power A^B (ln(A) db + B da / A)
 v Bth root BvA (da / A B - db ln(A) / B^2)
 L log base B (da / A ln B) - (ln A / ln B) (db / B ln B)
 # paste digits %%% reserved
 product rule d/dx f(x) g(x) = df g(x) + dg f(x)
 quotient rule d/dx f(x)/g(x) = (df g(x) - dg f(x)) / g(x)^2
 chain rule d/dx f(g(x)) = df/dg dg/dx
 An example of applying the above to derive the formula for the derivative
 of A^B (where A and B are both functions of x):
 A^B = exp(ln(A) B)
 d/dx exp(Y) = exp(Y) dY
 d/dx ln(A) = 1/A da
 d/dx ln(A) B = db ln(A) + B/A da (by "product rule")
 d/dx A^B = d/dx exp(ln(A) B)
 = exp(ln(A) B) d/dx (ln(A) B)
 = exp(ln(A) B) (db ln(A) + B/A da)
 = A^B (db ln(A) + da B / A)
 = A^B ln(A) db + B A^(B-1) da (standard form)
 d/dx ln(A) / ln(B) = (ln(B) da / A - ln(A) db / B) / ln(B)^2
 log(a+b) = log(a(1+b/a))
 = log(a) + log(1+b/a)
 Another example for atan2(a,b):
 d/dx(arctan(a/b)) = 1/((A/B)^2+1) d/dx(A/B)
 = 1/((A/B)^2+1) (B da - A db / B^2)
 = (B da - A db) / (B^2 * ((A/B)^2+1))
 = (B da - A db) / (A^2 + B^2)
 which agrees with Wikipedia (at en.wikipedia.org/wiki/Atan2#Derivative):
 Informally representing the function atan2 as the angle function
 theta(x,y) = atan2(y,x) [...] yields the following for the differential:
 dtheta = d/dx(atan2(y,x) dx) + d/dy(atan2(y,x) dy)
 = (-y/(x^2+y^2)) dx + (x/(x^2+y^2)) dy
 
It is fairly easy to show that in any equation that constitutes a
solution, if any expressions or subexpressons evaluate to 0, the
entire equation can be replaced with an equivalent form that does not
involve 0. The equivalent form is never more complex and is usually
simpler. It can also be shown that any solution involving a
subexpression that contains x and has a 0 derivative can be reduced
to a simpler solution that does not.
 Why RIES Calculates Derivatives
Consider the value X = 2.5063. Each of the following solutions (or an
algebraic equivalent) will appear when you run RIES on it:
 2 x = 5 for X = 2.5 *
 x^2 = 2 pi for X = 2.506628274631 *
 x^x = 1+9 for X = 2.5061841455888
 x^2+e = 9 for X = 2.5063356063267
The exact solutions to these equations are all different values,
of course, and they all form successively closer approximations to
the supplied value 2.5063.
%%% the following description doesn't match the numbers; don't know
how to fix it...
Look at the last two. Notice that the supplied value, 2.5063, is
between the exact solutions of these two. Also, the equations can both
be expressed in a different way:
 2 X = 5 instead of X = 5/2
 X^2 = 2 pi instead of X = sqrt(2 pi)
Consider what would have happened if the last digit had been one
higher or one smaller. How much does it change the "closeness" of the
match? With the original two equations, where there is just an "X" on
the left-hand side, this is easy to figure out: if X is 0.0001 lower,
it's 0.0001 closer fit for the first equation, and 0.0001 further from
the second. But look at the alternate forms: If you subtract 0.0001
from X, that subtracts 0.0002 from 2 X, but it subtracts 0.0005 from
X^2. That's a big difference -- if we're using these forms of the
equations (as RIES does) to determine how well the two sides match,
then this altered value of X makes the match "move" further with
respect to one equation than it does with respect to the other!
The reason this happens is, of course, because we're looking at an
expresson on the left-hand side, rather than just a single "X".
Putting an expression on the left-hand side makes it harder to see how
close a match you've got.
If you still don't believe this, consider 1.047246, and exclude sine
and cosine from the function set:
<pre>
 ries 1.047246 -NSC
 Your target value: T = 1.047246
 3 x = pi for X = T - 4.84488e-05
 x^5 = 3 root-of 2 for X = T + 4.81228e-05
 1/ln(sqrt(x)) = 4^e for X = T + 1.77179e-05
 ...
</pre>
We see that 1.047246 is an equally good solution for the following:
 x = pi / 3 (too high by 0.00005)
 x = 15th root of 2 (too low by 0.00005)
but if you express the solutions as the #ries# output does:
 3 x = pi (3.14172, too high by 0.00013)
 x^5 = cuberoot(2) (1.25992, too low by 0.00029)
suddenly it looks like the 3 x = pi solution is more than twice as
good. #ries# notices this and compensates for it regardless of the
form in which the equation is actually found. If you run RIES on the
value 1.047246 it will present both solutions, in the following form:
 X * 3 = pi
 X ^ 5 = 3 v 2
The philosophy adopted by RIES is that the "true" form for evaluating
the closeness of a match is the form where there is just one X, and
nothing else, on the left-hand-side of the equation, and just numbers
and symbols, but no X's, on the right-hand-side. (Let's call this the
"normalized form".) Put all equations into normalized form, calculate
the value on the right and look at the difference between X and this
value to determine how good the match is.
But now go back to our 2.5063 example above -- one of its solutions
was:
 X^X = 10
This equation *cannot* be reduced to something with just one X on the
left-hand side -- there is no "inverse-of-X-to-the-X" function fn[1]. What
does RIES do?
fnd[1] Actually, there is, but it uses an obscure function called the
"Lambert W function" which is the inverse of *y*=*xe^{x}*. More
details are [here|+numbers:xxy_root2] if you are interested.
RIES calculates derivatives. By using the value of 2.5063 for X and
calculating the derivative of X^X for this value of X, you get (about)
19. A derivative of 19 means that any small change in X will cause a
19-times-bigger change in X^X. That's important, because it allows us
to compare the closeness of the match on "equal footing" with other,
normalized equations like X = sqrt(2 pi).
Derivatives work so well, in fact, that RIES does not even have to
bother solving its equations for x. ven if it were easy to do this
(which it is not), leaving the equations unsolved is still an
important speed improvement. Part of the reason RIES is so fast is
that it generates left-hand-sides and right-hand-sides separately
(like a wl[Bidirectional_search] in graph theory) and tries all the
combinations to find possible solutions. It can do this a lot quicker
because a half-equation is smaller than a full equation, and therefore
there are less possibilities to check out.
Furthermore, derivatives allow RIES to quickly and easily check a
possible equation to discover the value that X would have to be in
order to both sides to match exactly (it is essentially performing one
step of Newton's method).
Examples of LHS and RHS expressions for the test case 2.5063. These
are shown in groups that correspond to the pairs that would actually
generate matches in a search:
 expression value deriv.
 5/2 2.5 0
 hyprt(10) 2.506184 0
 x 2.5063 1.0
 ,/25 5 0
 2 x 5.0126 2.0
 x^2 6.281540 5.01
 2 pi 6.283185 0
 x^2+e 8.999822 5.01
 9 9 0
 9+1 10 0
 x^x 10.00222 19.19
The derivative is based on the concept of an imagined error-bar in x
that is assumed to be small enough so that all reported matches are
relevant, but which is not zero. Thus, it is in units of the
infinitesimal quantity "epsilon".
Let us now imagine that there is a constant "g" equal to (pi-1) *
2.5063 ~= 5.367473. Then we would have the following grouping:
 expression value deriv.
 pi^2-2 7.869604 0
 x+g 7.873774 1.0
 pi x 7.873774 3.14
The matching "pi x = pi^2-2" constitutes a closer match than "x+g =
pi^2/2" because in the former case, x would only need to be decreased
by about 1/pi as much to make it an exact match. So, the closeness of
a match is measured as |LHS-RHS|/derivative, where smaller is better.
To enable actual error bars to be provided with data, a value of
"epsilon" must be adopted for use with notionally precise supplied
values. This can be gleaned from the number of supplied digits in the
input, or the precision of the floating point format can be used. In
the latter case we would take the data value divided by 2 to the power
of the number of digits in the mantissa.
 
Future enhancements:
 See "UNFINISHED WORK" section above
 
REVISION HISTORY
 20000207 Begin (it does not do much except parse the parameters)
 20000208 Add parsing of level adjust and more design comments
 20000209 Write code to generate forms. "Discover" the Motzkin numbers.
 20000210 Optimize gf.1() by making it compute the stack depth as it
goes along rather than repeating the whole stack history on each test.
This improves time to generate all forms of length <=19 from 188
seconds down to to 37.9 seconds. Add tracking of min and max weights,
and start writing expression generator.
 20000215 It now generates forms within min and max possible
complexity limits, to avoid generating expressions from forms that
cannot possibly fall within the current complexity limit. Add a little
optimization; speeds up a deep search from 53 seconds to 37 seconds.
 20000217 It generates expressions, but doesn't prune for complexity
or evaluate. For the initial groups of 5, the numbers of expressions
evaluted are: {0, 12, 12, 72, 1368, 1368, 17928, 345672, 5875272,...}
 20000217 Prune for complexity limits; each expression is now
generated exactly once. The numbers are down to: {0, 1, 3, 18, 69,
182, 1046, 5358, 27123,...} Still need to prune foolishness like "11+"
and "nn".
 20000217 Prune almost all obvious trivial patterns. The numbers are
now down to {0,1,3,13,43,122,486,2186,9775,...}. Then prune [JK+] and
[JK-] for small integers, [jK*] and [jK+] for any j<K, and a couple
lesser things, and they're down to {...,43,106,391,1861,8608,...}.
20000217 Write metastack routines and exec(), and prune on eval errors
like divide-by-zero. This brings the numbers down to
{0,1,3,13,41,96,336,1605,7129,...}
20000217 Add using.x (formerly "on_rhs") variable and AM_RHS so it
generates LHS and RHS expressions, and make main outer loop call
gen_forms twice. In the LHS case there is lots of pruning. Figure out
that there is a substantial problem with error factors (derivative
of X) that will make lots of bogus "solutions" show up in the output.
20000218 Finish figuring out how to deal with the bogus solutions --
compute derivatives on all LHS expressions. This also solves the
"trivial solutions" problem (e.g. "x - x = ln(1)"
20000220 Add derivative calculation to exec() and the ms_xxx
routines. Initial output looks good, but I'm a bit worried about the
very high derivative values on things like [x8s^] -- they might
cause lots of bogus matches on things that aren't even close, like
"x8s^ = 6".
20000220 Add pruning of zero subexpressions, and of LHS expressions
with zero derivative. Examine output to check this; discover a bug
that ultimately turns out to be because I had ">> 2" instead of ">> 1"
in the s[] and ds[] declarations in struct metastack.
20000220 Refine the termination condition (now it counts generated
expressions, rather than cutting off at a certain complexity score --
this is to gain independence from the specifics of the weights) and
discover another dimensioning bug in struct metastack. Benchmark
changes in PASS_GRAN: When it is set to 4, 2, and 1 the execution time
for generating all expressions up to complexity 64 (2236462
expressions) is 3.57, 4.16, and 5.24 respectively. This actually a lot
better than I thought it would be -- I thought there would be a lot
more overhead from repeating the same subexpression evaluations over
and over again.
20000220 Make it increase LHS and RHS complexity limits at independent
rates, such that the population remains equal between the two sides.
Add a bunch of comments documenting what's been done so far.
20000221 Write bt.insert(); it now reports exact matches! It finds,
for example, [x3+] = [4s] for X=13. However, for X=143 it reports
[xrx*] = [1], with a derivative of 1.9e-19. Add PRUNE_DERIV to try to
fix this, and it starts reporting [xxqL] = [2] with a derivative of
0.00258248 -- this turns out to be a bug in the deriv formula for 'q'.
Eventually increase PRUNE_DERIV to around 1e-14, then decide to make
it a variable and add p_ovr and n_ovr. Write bt_prev() and bt_next(),
but not using them yet.
20000221 Figure out that I can check just the nearest LHS on either
side of an RHS, and vice versa. Write check.sides() and check.match.
It now finds answers and prints them out!
20000221 Fix bug in exact-match reporting. When given 1.5065916, it
reports:
 match: [xe+s1+] = [p6*] (solution is X+5.14855e-08)
Make it report delta you'd have to add to X to get each match to be
exact. Test cases that currently produce bogus results: 'ries -l1 27',
'ries -0.28676844', 'ries 403'
20000222 Adjust best.match by 0.999 each time to avoid long strings of
roundoff-error results. This fixes the './ries 403' case. Rename to
"ries" (it used to be called "misc"). Start analyzing memory usage.
Figure out how to save some memory in the node structures, and more
importantly, how to group nodes together in physical memory such that
the program degrades more gracefully when physical memory limits are
exceeded.
20000222 Implement -S and -O command-line options. This has the
side-effect of making certain bogus match bugs easier and quicker to
reproduce.
20000222 Change -S and -O to -N and -S respectively.
20000222 Convert manpage to #nroff# format. Improve implementation of
-S/-O/-N precedence. Implement -O option. Figure out why restricted
symbolsets create very small numbers of equations (the level controls
the number of expressions generated, which is always larger than the
number of expressions inserted in the tree)
20000222 Fix bug that caused negative X's to give no solutions -- it
was initializing best.match to a fraction of X, and not taking the
absolute value! Implement -i option so I can find an expression for
70458.
20000223 Ignore -i if target is non-integer. Adjust weights for seft
'b' operators; 'l' now appears much less often in expressions. Add -y
option.
20000223 Write add.rule(); convert all AM_xx rules to add.rule()
calls. Now the AM_xx rules degrade gracefully with the symbolset
options. Change -y option to -x.
20000223 Clean up formatting in report.match(); add complexity score to
output. Write perl script to benchmark and take statistics on -i option
for a wide rance of integers; leave it running overnight. Change all 'int'
to 's16'.
20000224 Today's date as a mathematical expression:
 (((4^(4^(1/e))-pi)^2)-pi)^2 the repetition is cool.
20000224 Write initial version of infix.1, then add a few rules
(ordering of bare symbols in '+' and '*'; always reverse order for 'v'
and 's'). Add -F option.
20000224 Add memory usage statistic.
20000225 Write infix.preproc; implement parentheses precedence.
20000225 Write eval() and newton(); printed X values are now (almost
always) exact roots of their equations. Add AM_1K rules.
20000226 Add copyright and GPL notices; add URL to first printf
20000226 Add symbol definition strings
20000227 Add sin and cos operators. Eliminate several bogus exact
match errors related to roundoff and loss of accuracy, e.g.
cosine(0.0001). Fix major bug in k_prune.deriv test: it only pruned
positive small derivs, not negative small derivs. Add symbol FORTH
names and write postfix.formatter (but it isn't used yet). Don't report
matches if newton() returned an error.
20000227 Add -ie variant to -i option (only_exact)
20000228 Kill another missing-fabs bug; remove prune on "almost exact
matches" which is now adequately covered by the derivative tests.
20000228 Add 'E' operator, AM_l and AM_E attributes, and PS_REVPOW.
check.match now uses newton() to evaluate score more accurately.
Add loss-of-significance test to '-' operator.
20000228 check.match calling newton was royally slowing things down.
Now it uses the old, much quicker test and then uses newton() as a
confirmation test. This cuts time for "ries 2.5063" from 4.09 down to
1.39 on the Cyrix 180 (which is about what it has been since 0223).
Add statistics of pruned subexpressions ("dead-ends").
20000229 Fix bug in newton() that prevented success if target was
negative. Add debug printf's 'mnr'.
20000229 Add time display; add debug printf's 'opqsABCDEFG'.
20000229 Add debug printf's 'Hy' and a few notes about complex analytic
definitions
20000301 Change debug printf 's'; add another loss-of-significance test
to '+' and '-'; add debug printf's 'IJKLtuvx'. Add 'I' operator and
special-case tests for no defined symbols of each seft (this is a massive
optimization for "ries -ie 7 '-S1+*-/^v'")
20000302 Add debug printf 'w' and improve 'r'. Add arctan function (but
not using it yet) and some notes about derivatives.
200003xx Benchmarks on a 333-MHz Celeron:
 -command-------------------------- -mem- time equations
 ries 2.5063141592653589 960K 0.3 9.5619e7
 ries -l1 2.5063141592653589 3072K 1.5 1.0306e9
 ries -l2 2.5063141592653589 10.7M 6.5 1.2832e10
 ries -l3 2.5063141592653589 32.7M 24.4 1.2042e11
200203xx Benchmarks on an 800-MHz PowerPC G4 (iMac, model M6498LL/A):
 -command-------------------------- -mem- time equations
 ries 2.5063141592653589 960K 0.1 9.5062e7
 ries -l1 2.5063141592653589 3072K 0.8 1.0248e9
 ries -l2 2.5063141592653589 10.6M 3.6 1.2765e10
 ries -l3 2.5063141592653589 32.6M 14.3 1.1976e11
 ries -l4 2.5063141592653589 114M 70.3 1.4609e12
20020610 Figures from an unknown test (I can't find a record of the
details). I think the numbers are: memory usage; total
expressions/total distinct; total equations.
 -l0 1136K 32198/14576 53165000
 -l1 2876K 145183/51492 666660000
 -l2 7632K 530668/152904 5712000000
 -l3 27100K 2311204/559419 78120000000
200310xx Benchmarks on an 800-MHz PowerPC G4 (iBook G4, model M9164LL/A):
 -command-------------------------- -mem- time equations
 ries 2.5063141592653589 960K 0.1 9.5062e7
 ries -l1 2.5063141592653589 3072K 0.7 1.0248e9
 ries -l2 2.5063141592653589 10.6M 3.2 1.2765e10
 ries -l3 2.5063141592653589 32.6M 12.8 1.1976e11
 ries -l4 2.5063141592653589 114M 62.3 1.4609e12
 ries -l5 2.5063141592653589 398M 388 1.7825e13
20050715 Benchmarks on a 2-GHz PowerPC G5 (part of a dual system,
model M9455LL/A):
 -command-------------------------- -mem- time equations
 ries 2.5063141592653589 960K 0.0 9.5062e7
 ries -l1 2.5063141592653589 3072K 0.2 1.0248e9
 ries -l2 2.5063141592653589 10.6M 1.5 1.2765e10
 ries -l3 2.5063141592653589 32.6M 6.6 1.1976e11
 ries -l4 2.5063141592653589 114M 36.2 1.4609e12
 ries -l5 2.5063141592653589 398M 260 1.7825e13
 ries -l6 2.5063141592653589 1.41G 2183 2.2297e14
20070511 Alan Eliasen emails me to tell me that 'ries -l4 193707721'
gives the spurious result "x.(1/S(p))/x = 1/S(p) (exact match) {107}".
This is a problem I have known about for a while. I reproduce the
problem and reduce the symbol set to make it appear faster; it can be
reproduced with:
 ries -F '-SxpSr/ *' -l0 1002353667
Add rule ("",'S', AM_pi, 0), which eliminates "sin(pi)". It now reports
the following more complex version of the same bug:
 x.(1/S(p.p/p))/x = 1/S(p.p/p) (exact match) {181}
{This isn't actually a problem, it's doing what it should: When the
symbol '1' is not allowed, then "K/K" for any constant K *is* allowed.
However "ries -F '-SxpSr/ *1' -l1 1002353667" does produce the bug.
See 20090513 for fix. -20090513}
20090216 Clean up the comments, and add the "ARCHITECTURE" section,
plus a few notes on how to approach a multithreaded implementation.
 Here are the benchmarks from the 733 MHz Pentium 3 (from the manpage,
which I have now decided to update):
 memory equations digits runtime
 usage tested matched (733MHz P3)
 -l0 960K 95,000,000 6+ 0.1 sec
 -l1 3.1M 1,030,000,000 7+ 0.7 sec
 -l2 11 M 12,800,000,000 8+ 3.3 sec
 -l3 33 M 120,000,000,000 9+ 12 sec
 -l4 115M 1,470,000,000,000 11+ 63 sec
and here are the benchmarks from the 2.33-GHz Core 2 Duo (MacBook
Pro):
 -command-------------------------- -mem- time equations
 ries 2.5063141592653589 960K 0.0 9.4018e7
 ries -l1 2.5063141592653589 3008K 0.1 1.0153e9
 ries -l2 2.5063141592653589 10.6M 0.6 1.2643e10
 ries -l3 2.5063141592653589 32.4M 2.7 1.1865e11
 ries -l4 2.5063141592653589 113M 13.6 1.4461e12
 ries -l5 2.5063141592653589 396M 81.7 1.7642e13
 ries -l6 2.5063141592653589 1.40G 651 2.2062e14
20090510 Add standard output format, and make it the default.
"-F" now takes a numeric argument (which formerly it did not)
and "-F" alone defaults to "-F0" which is what it did before.
{Formerly "F" was for "FORTH", and now it stands for "Format".}
 Future idea: A -W option to view weights (same as -S) or define
weights individually. --include could then be used to select weight
presets (e.g. a set of weights for electrical engineering). {This was
implemented on 20121209}
20090511 In infix modes, display operator 'L' (log base A of B)
in advance of both of its arguments. No longer put parens around
single-symbol argument of negate (e.g. emit "-x" instead of "-(x)").
20090513 Add rules to put bare 'x' after non-x-containing expressions
and bare 'pi', 'phi' and 'e'. Finish postfix.formatter (which had
never been brought into spec with the other formatters) and add it
as a 4th output option.
 As noted at 20070511, there has been a problem with expressions like
"pi*pi/pi". The command "ries '-SxpSr/ *1' -l1 1002353667 -F0" currently
shows the problem. I fix this by adding the three-symbol pattern AM_KxK
and the rule ("", '/', AM_KxK, 0).
20090515 "ries -l6 2.5063141592653589" running uncontested on the
2.26-GHz Nehalem uses 1.40G of memory and tests 2.2154e14 equations in
442.4 seconds. Another test: 1.3525746932102463: 1.43G, 2.3011e14,
463.2. These times are 47% and 40% faster than the Core 2 Duo.
 Slight improvements to -Ds output. Improve rules for emitting "*",
" ", or "" for multiplication in different situations. Add rule
("12345678n", '-', AM_jK), and def_amkey/sym_amkey[] optimization.
20090808 Remove "val" parameter of my_alloc()
20101218 Increase default level to -l2 (while preserving the effect
of all '-l' values when such an option is provided: "ries 1.234" is now
the same as "ries -l2 1.234", but "ries -l1 1.234" does the same thing
it used to do).
 Print 'x' instead of 'X'; use strcmp and strncmp in a few places.
20111216 Increase MAX_ELEN from 16 to 19 since we have room (struct
size is 64 in either case). Comment out unused obt_xxx declarations.
20111218 check.sides now checks more than one neighboring LHS if the
newly-inserted node is an RHS. Oddly, this does not cause more matches
to be reported, and I am not sure why.
20111219 print.end now displays max LHS and RHS complexity values; add
--find.expression option.
 Add k.vanished_dx and use it in several places to enforce a stricter
constraint on derivatives. This eliminates many of the more obscure
tautologies like "x^(4/ln(sqrt(x)))".
 Add --eval-expression option, and error numbers and error strings to
support it.
20111220 Restore ERR_EXEC_TRIG_LOW_DX error for RHS expressions;
restore legend in normal infix display format; eval() now reports
ERR_EVAL_TOO_LONG.
20111221 Slight improvements to postfix.formatter; report (most)
unrecognized or badly-formatted command options.
 Catch a few more tautologies by comparing magnitude of dx to
magnitude of x.
20111222 Add --version option.
20111223 -S without any symbols displays the symbol table and exits.
20111224 Enhance -Ds output by explicitly mentioning the Newton-Raphson
step.
20111226 Add loss-of-significance tests to each of the operators that
seem to need it: ln(x), e^x, sin, cos, +, -, a^b, a,/b, and log_a(b).
 Significance loss errors have always been easy to find in ries
output, by giving a target value that has a simple solution (like
"ries -1.4142135", which is used throughout this note). After ries
gives the simple solution (in this case "x^2 = 2"), subsequent
solutions can incorporate the simple solution in a
loss-of-significance tautology.
 The first of these solutions was {2012.0505: Now considered okay, now
that addition actually checks for loss of precision explicitly}
"1/(x+sqrt(2)) = e^(4^2) {86}". This incorporates the direct solution
into a tautology by turning "x^2=2" into "x+sqrt(2)=0", which isn't
exactly 0 because x is not exactly sqrt(2). Significance can be said
to be lost because many of the significant bits of the two terms "x"
and "sqrt(2)" cancel each other out in the addition. {2012.0505: In
this case, handled in exec() case '+', the loss of significance in
addition is easy to test rigorously: if (a+b)-b is equal to a, and
(a+b)-a is equal to b, then no information has in fact been lost.}
 This "x + -x" case is probably the most common source of significance
loss that I had not addressed until today. I fixed it by adding the
"fabs(rv) < (fabs(a) * k_sig_loss)" tests in exec() label
"add_common:"
 "ries -1.4142135" then gave "-ln(x^2-1) = 1/7^8 {109}" until I added
the k_sig_loss test in exec() case 'l'
 It then gave "-sin(pi x^2) = (1/5)^9 {111}" until I added the
k_sig_loss test in exec() case 'S'
 It then gave "log_2(x^2-1) = 1/-(e^(e^e)) {117}" until I added the
k_sig_loss test in exec() case 'L'
 Then it gave "1/(1/(x-1)-x) = (e^(4^2)),/2 {122}" until I added the
k_sig_loss test in exec() case 'v'
 It then gave "(x-1/(x+1))^2 = e^(1/(e^7)^2) {125}" until I added the
k_sig_loss test in exec() case 'E'
 After doing all of the above there was only one operator ('^') that
seemed to need a k_sig_loss test.
20111227 Change many details of formatting and printing to accomodate
the recent changes in handling of significant digits. Add constants
k_nominal_digits and k_usable_digits and associated format strings
fmt_g_xxx and use these in most places a floating-point value is
displayed. Redo the equation justify (space-padding) code to make the
best use of 80 columns while still showing centered equations and 15
significant digits when the -x option is given.
 Add init.numerics(); k_phi and k_pi are now computed from scratch.
 New 'z' option to -D and debug_z flag to show messages printed by
init.numerics().
20111228
 Add more k_sig_loss tests in sin and cos.
 Display CPU time as "%.3f" rather than the old "%d.%d" that showed
seconds and tenths.
 RIES now exits when (best.match < k_0) regardless of the got.exact
flag; this fixes a bug that would cause RIES to loop forever if it had
not yet gotten an "exact match" at the point when best.match goes
negative.
 init.numerics() now computes k_e and detects the size of the ULP
(Unit in Last Place or "least significant bit").
 Add --min-match.distance option and g_min.matchsize to prevent
reporting of really close matches (useful in e.g. finding classical
approximation formulas for pi)
 Add --significance-loss-margin and k_sig_loss to allow going back
to RIES behaviour before all the new k_sig_loss tests were added.
 Share code by combining two exit tests into exact_exit() routine.
 ln(a) and log_a(b) functions now no longer reject a case like
ln(1.23e10) (in which only about 1 significant digit is lost)
 k_prune.deriv had been used for three different purposes, and now is
being replaced with three separate variables. The three purposes of
k_prune.deriv were:
 1. Detecting convergence in newton(). Now using new k_newton_settled
 2. Pruning subexpression tautologies in ge.2(). Now using k.vanished_dx
 3. Setting n_ovr and p_ovr in init2(). For now this is unchanged.
20111229 The big trig change: S and C now take units of pi radians, so
"3rS" now produces the value sin(pi/3)=0.866025... To support this I
also add the --trig-argument-scale option, which if given with the
parameter "1" gives the old trig units (radians). As a bonus, users
can now get degrees or grads or any other angle units.
 To avoid confusion, the trig functions are now called "sinpi" and
"cospi" in the output, unless --trig-argument-scale is set to
something other than pi; and the function definitions make this
obvious by adding a message like "For the trig functions, 360 units
are a full circle." at the end.
 While I'm at it, I add the tangent function 'T'. This had been in
the manual and in the symbolset since the start, but had not yet been
written in exec(). This of course makes a lot of things solve more
quickly, such as Gosper's "0.9674026381747" (the root of "tan(x) =
3x/2", which previously needed a -l5 search, but with the unit change
and the tangent function, a -l0 search is sufficient.)
 One result of the trig units change is that there are a lot more
natural identities connecting common fractions and small radicals. For
example sqrt(2) is 2 sinpi(1/4), sqrt(3) is 2 sinpi(1/3), etc. This
makes ries generate and insert fewer expressions at a given complexity
level. It's so good in fact that you get higher complexity *and* less
memory usage even after adding the tangent function. Here are some
benchmarks, showing how much memory it used up and what level of
complexity it was able to achieve:
 --------old-------- --------new--------
 complexity complexity
 target-value LHS:RHS memory time LHS:RHS memory time
 1.5063 -l2 66:61 13632 KiB .319 66:62 12928 .342
 .328106566874978253 80:76 486848 KiB 39.5 81:76 461312 38.0
 1.9511889024071 80:76 494016 KiB 39.7 81:76 458432 37.8
 .922524879060934752 80:76 499264 KiB 40.4 80:77 464320 40.8
(All except the first is at level -l5; the old figures were measured
using the options "--trig-argument-scale 1 -NT")
20111230 Add a pruning rule for [K+K-]->[] and pattern AM_KpK (which I noticed
after updating the sin(4)/cos(4)=1.1578212823 example from the manpage
and discovering that "ries 1.1578212823" prints "x+1-1 = tanpi(4/pi)".)
 Time measurement is now done by wall-clock time rather than CPU
cycle time. We do this with the gettimeofday() function. Formerly we
used getrusage(), which under-reports elapsed time on Nehalem and
later Intel microarchitectures when the core that we are running in
also has another thread running. This change also involved adding the
inittime() routine.
 Make several changes to facilitate porting this file to other systems:
 * #include <limits.h> and use LONG_MAX instead of __LONG_MAX__ (which
was a GCC-only predefined constant)
 * The symbol stack effect attribute (formerly called "class") is now
called "seft" (both in the comments and more importantly in variable
names) to facilitate porting to C++ (where "class" cannot be used as a
variable name because it is a language reserved keyword)
 * As mentioned above we now measure time with gettimeofday(); in
Windows we provide a gettimeofday() based on _ftime().
 * #include "stdafx.h" has been added before the other includes, to
facilitate Visual C++
 * #define RIES_VERSION ... has been moved to fall after the includes.
20111231 Add bt_first(), not yet used.
20120102 Increase k_min_best.match when fabs(target)>1.0, to avoid a
"search forever without exiting" condition that would always happen when
fabs(target)>16.0
 Add --include/-p (load profile) option, implemented in new routines:
file_read, delimit_args, and parse.args (which includes all the
argument-parsing formerly in main()).
 An included file contains options, delimited by blank space with
optional comments starting with '#'. Each non-blank non-comment
'token' is treated as if it were a string in argv[], with the strings
inserted into argv[] at the point where the settings-file option is
given, with recursion to support nested includes. Max 10MiB per file,
max recursion depth is 27 levels. This is to provide for all the fancy
new features I'd like to add someday (like user-defined symbol-weights
and functions).
 Add bt_depth and bt_stats (which is run by debug_y) to see if
rebalancing the binary tree will help. I conclude that it will not
help much: in a typical run the average tree depth with rebalancing
would be 17.5, and the actual (unbalanced) average depth is 24.5.
20120103 Add check.exact_match(), removing more RIES-specific code
from bt.insert.
 Add next_isparam(), clean up argument parsing and fix a bug with
parsing --significance-loss-margin argument.
20120104 Improve several of the debug_n printfs; add a bunch of
debug_q printfs to answer the question of why the check.sides change
of 20111218 did not produce more matches.
 Rename the old "on_rhs" variable to "using.x".
 Remove a lot of obsolete comments relating to an old idea for
maintaining two separate (LHS and RHS) lists. Update the comments on
derivative formulas and the Gamma function.
 check.sides now always checks only one expression to either side of
the newly-inserted expression. This causes a very few changes in
output, mainly in the first one or two reported matches. The reasons
are explained below in a block comment "THE CHECK_SIDES PARADOX".
 To support the check.sides investigation I added the -D0 option and
debug_0, which does a complete dump of the expressions database after
every gen_forms pass.
20120105 Add stack overflow and underflow checks in eval(). Increase
MAX_ELEN again to 21 (the expr struct is still 64 bytes).
 Add MAX_SEFT_POP checks in add.symbol().
 Add better documentation of the workings of the metastack routines.
 Implement canon.val() but leave it disabled by default right now.
This tries to transform expressions into forms that have a value in
[1.0,2.0), which also transforms expressions. For example, compare:
 ries 1.50631415926535897932 --canon-reduction 0
 x-1 = 1/2 for x = T - 0.00631416 {50}
 1/ln(x) = sqrt(6) for x = T - 0.00213357 {62}
 x^4 = 2+pi for x = T - 0.000489459 {68}
 (x+1)^2 = 2 pi for x = T + 0.000314115 {69}
 x^pi = tanpi(sqrt(2)) for x = T - 0.000224565 {70}
 x^2 = tanpi(1/e) for x = T + 9.42083e-05 {60}
 x^2-2 = 1/(1+e) for x = T - 1.35846e-05 {79}
 1/(log_3(x)) = 3-1/pi for x = T + 8.87321e-06 {87}
 1/x+sqrt(2) = 1/ln(phi) for x = T - 2.77223e-06 {86}
 1/(log_7(x)) = 5-1/4 for x = T - 6.29902e-07 {94}
 . . . . . . . . .
 log_(2/sqrt(3))(x) = 3-1/2^e for x = T - 1.21748e-09 {128}
 max complexity: 66 62 128
 dead-ends: 2290352 5560466 7850818 CPU time: 0.369
 expressions: 183663 420208 603871
 distinct: 92770 113599 206369 Memory: 12928KiB
 ries 1.50631415926535897932 --canon-reduction nr2
 x = 3/2 for x = T - 0.00631416 {48}
 (1/ln(x))/2 = sqrt(6)/2 for x = T - 0.00213357 {98}
 x^4/2 = (2+pi)/2 for x = T - 0.000489459 {104}
 (x+1)^2/2 = 2 pi/2 for x = T + 0.000314115 {105}
 -(x-4)/2 = pi,/2 for x = T - 5.21371e-05 {95}
 (1/(x^2-2))/2 = -(pi-5) for x = T + 2.16534e-05 {110}
 (1/(x^2-2))/2 = (1+e)/2 for x = T - 1.35846e-05 {115}
 e^(x-1) = (1/pi+3)/2 for x = T - 5.75694e-06 {101}
 1/(x/phi)^2 = 8,/pi for x = T + 7.44845e-07 {94}
 sqrt(3-x) = -tanpi(e-1) for x = T - 6.44308e-07 {101}
 (1/-(log_7(1/x)))/2 = (5-1/4)/2 for x = T - 6.29902e-07 {144}
 . . . . . . . . .
 sinpi(1/5),/x/2 = 1/-cospi(1/ln(sqrt(7))) for x = T + 1.15156e-10 {151}
 max complexity: 66 62 128
 dead-ends: 2290352 5560466 7850818 CPU time: 0.369
 expressions: 183663 420208 603871
 distinct: 63206 67234 130440 Memory: 8192KiB
The intent is to get further with a given amount of memory by
collapsing equivalent equations like "x=-2" and "-x=2" into the same
solution. It succeeds in doing this, but also causes different matches
to be discovered and reported (despite the fact that the main loop is
still making the same choices about whether to call gen_forms on LHS
or on RHS). I need to investigate this further to discover why (one
possibility is that the new solutions come from things that would have
matched earlier, but didn't, because they lie in different places on
the number line and wouldn't meet up without needing more complexity).
20120107 -pfoo now tries opening "foo.ries" if "foo" cannot be opened.
 Add thrash_check() and related real-time memory allocation
benchmarking (not yet finished). This is based on ideas from a
discussion with Charles Greathouse yesterday.
20120108 More work on thrash_check(); add --min-memory option (which
presently does nothing, but eventually will prevent RIES from quitting
early in the event of unanticipated system slowness). Add --max.memory
option (which makes RIES definitely quit before exceeding a given
amount of memory).
20120113 Move nested functions out of parse.args() for compatibility
with non-GCC compilers.
 Presently, "ries 0.9674026381747 --eval-expression xp/T" shows a
different answer on different computers: I get "tanpi(x/pi) =
1.4511039572620501" on MP16 and "tanpi(x/pi) = 1.4511039572620503" on
MBP. This causes commands, such as "ries 0.9674026381747 -l3", to
produce different results ("x/tanpi(x/pi) = 2/3" on MP16 and the
bizarre "x/tanpi(x/pi) = log_(sqrt(8))(2)" on MBP).
 I suspect a variation in the maths library routines under different
versions of MacOS. To fix this I probably need to bring in my own
tangent function. Sources for one version are in
"sun-trig-functions.txt". {This was addressed on 20121215, when I
added msal_math64.c and the RIES_USE_SA_M64 switch}
20120115 Add --memory-abort-threshold option and fix some segfaults
that were happening if the user does not give an expected option
argument.
 Add typedef sym_attr_block and use it for all the symbol table arrays.
20120121 Move all symbol tables into a single array of sym_attr_block
20120122 Add CANONVAL_MUL2 operation and stub cv.simplify() routine.
20120423 Notes on what I need to accomplish to finish canon.val,
decanon, and cv.simplify routines:
 canon.val adds 'n', 'r', and/or '2/', '4/' to try to put an
 expression's value into the range [1, 2). Call it just before
 calling bt.insert. It should work on the existing expression
 structure and metastack state, and it needs to stop if MAX_ELEN is
 reached. We also need a decanon routine to undo its work before
 proceeding with the rest of ge.2.
 cv.simplify operates on two expressions (an equation), removing
 any unnecessary symbols that may have been added by canon.val. It
 takes two expressions A and B and removes any trailing '2/', '4/',
 'r' and 'n', or possibly shifts symbols from one expression to the
 other, all so as to make the equation A=B lower-complexity. Call
 this in check.match right after the initial test passes and before
 running newton() on the equation. This will most probably involve
 having two dedicated expression data structures so it (and newton(),
 etc.) will see our simplified equation rather than the two
 expressions that are actually in the database.
 The purpose of this is to conflate values like -2, -1, 1/2, 1, 2
 into a single tree entry, which should make more effective use of
 memory, at the expense of the output being a little less
 well-ordered by complexity-score.
 {These changes were begun on 20120514, but not finished until bug
 fixes on 20121210.}
20120425 RIES was used by Randall Munroe for XKCD #1047 (xkcd.com/1047)
{This led to changes allowing RIES to function as the back-end to an
online RIES server.}
20120428 Changes for UNIX and Windows compatibility (suggested by
Markus Milleder, via the xkcd forums).
 Also change Nth root symbol ',/' to '"/' because the '"' looks more
like the superscript 'n' that I use in the HTML for the ORIES server.
The server is mostly functional, at mrob.com/ries
20120503 Add --wide-output option (mainly for use by the ORIES
server). This is a 132-column version of the standard 80-column
output, adding two new columns. It displays the normal output plus the
-x version of the column showing equation roots, and a new column
showing the ratio (target/delta) in the form "1 part in 50" as in
xkcd.com/1047
20120505 There has been a bug for a while, causing "ries -l-2 -i 143"
to find the bizarre equation "(x-2)+3 = (3*4)^2" rather than the much
more obvious (and lower-complexity) "x+1 = (3*4)^2". This resulted
from sig-loss pruning and the k_sig_loss value 0.01: Because 1 is less
than 143*k_sig_loss, RIES disallowed the addition x+1. Similarly,
"ries 7775" failed to find "x+1=6^5" or even an alternate like
"(x+9)-8=6^5" because it couldn't add 1 (or any integer) to an x that
large.
 As a "narrow, conservative" fix, I have altered exec() case '+' to
attempt "(x+1)-1" and see if it is precisely equal to x. If so, then
we can assert that no significance loss has actually occurred
 Also add parse_target() to support future --mad option.
20120508 Numerous small changes to avoid warnings about sign mismatch
and using 'char' as an array index. Also add functions sym.strsym and
sym.strcmp.
20120509 Add --max.match-distance option, bringing it a bit closer to
supporting --mad.
20120510 Add the symbol ' ' which does nothing, so that --eval-expression
expressions can be formatted with whitespace if desired.
 Display a "still looking..." message if no results have been reported
after the first 2 seconds.
20120511 Add --numeric.anagram option and get it basically working. It
still needs a lot of testing, and the loop exit condition based on
search level needs a lot of work.
20120514 A bunch of little changes to fix compiler warnings.
20120515 Handle the gen_total / searchmax limits differently for
unidirectional searches like that done by --numeric.anagram. Prior to
this change, "ries -i 143 --numeric.anagram 143" would search forever
and finds nothing.
20120516 Add pf_float_wid; prune_count etc. are now doubles. Main loop
termination tests are now more flexible and general-purpose (important
for --numeric.anagram and other restricted searches).
20120518 Eliminate use of "s64" (long long 64-bit int) and associated
header files; I can use double instead for my large integers. Clean up
some size_t and char * typecasting. This is all for compatibility with
non-GCC compilers.
20120520 More cleanup, thanks to suggestions by Markus Milleder.
20120522 Add boolean datatype; eliminate all uses of explicit 32-bit
integers; add validate_types() and smarter detection/definition of s16.
20120613 Add --derivative.margin option. Also test different values of
the initial k.vanished_dx setting, and discover that it's safe to use
a much lower value like 1e-10. It's possible that most of the problems
fixed on 20111219 were also fixed by the sig_loss tests added later
in 201112; I'll have to do a lot of testing to be sure.
20120720 Add --one-sided option; --numeric.anagram now implies
--one-sided.
20120725 Change first argument of sym.strsym() from 'expr' to 'exp1' to
avoid namespace conflict with typedef 'expr'
20121202 Better implementation of "--max.match-distance 0"
 Implement the --match.all-digits option, which makes RIES work more
like ISC. The capabilities has already been implemented, primarily in
parse.target(), the --max.match-distance option, and the "Still
searching" message in the main loop.
20121205 If the target value is large in magnitude, automatically set
k_vanished.dx as if the --derivative.margin option was given.
20121206 Add several formatting routines (symstrncpy0, sym.strncat,
endstack, expr_break, expr_print_infix, eqn_print_infix) and
try.solve() to implement the --try-solve-for-x option.
 This is in a rather rough state, as exemplified by the command "ries
0.8183431428522 --try-solve-for-x", which will happily report "x =
sqrt(phi+pi)+3". I can fix that within try.solve() by calling eval()
after taking the square root of both sides to see if the LHS
is negative, and if so, append [n] to the RHS.
20121207 Fix the square root sign problem with 0.8183431428522, and
add some notes to try.solve.
 Increase NEWTON_MAX_ITER, and add notes about chaotic oscillation:
 2012.1207: "ries 0.2322795271987 -l0" does not find "sinpi(x) = 2/3"
 because Newton does not converge: sin(x) in that region causes a
 divergent oscillation. "ries 0.23227952719876987" works because it
 happens to hit on an exact match, which thus does attempt to call
 newton(). The oscillation can be seen by passing -Dn: RIES is clearly
 trying all sorts of equations involving trig functions like [xS],
 [x2-T], etc. I could "fix" this by making the Newton iteration only go
 halfway to the next point each time:
 curr = curr + 0.5 * ((rhs_val - lhs_val) / (lhs_dx - rhs_dx));
 and increasing NEWTON_MAX_ITER to about 100, but that also introduces
 bogus solutions like "tanpi(-x) = 7^2" that I might not want. It's
 clear that trig functions in x are being rejected fairly often and if
 I change this behaviour, it will change a lot of RIES output.
20121208 Fix derivative formula for sine (case 'S' in exec()). It
has been wrong ever since I added --trig-argument-scale, and was
causing newton() to not converge in many (if not most) cases.
20121209 Add period-2 loop detection in newton(); this might allow
eliminating k_newton_settled later.
 Add sym.strlen, sym.strtrail, bothtrail and sym.strclip and begin
implementing cv.simplify.
 debug_p replaces debug_q for the "first score not good enough" case.
 new flag debug_e replaces debug_E; debug_E replaces debug_F.
 debug_F replaces canon.val cases of debug_G
 add new flag debug_Q for cv.simplify
 make debug_S and debug_s distinct variables (but 'S' and 's' debug
functions still print the same messages)
 Fix derivative formula for cosine.
 report.match now takes symstr arguments; copy expressions before
calling cv.simplify (which fixes heisenbugs that were caused by it
modifying the expression strings in the actual database nodes).
 Add a couple more rules to infix.preproc to improve --canon-reduction
output.
 new flag debug_N replaces debug_o.
 debug_o replaces debug_p (including the debug_q case changed earlier today).
 debug_p is now used by infix.preproc
 debug_F and debug_f are now distinct (canon.val for LHS and RHS respectively)
 Fix derivative formula for tangent and increase its symbol weight to 6
(sine and cosine are both 3)
 Add --symbol.weights option. A few simple tests like "ries 2.5063
--symbol-weights 1:T" quickly reveal that these "Gosper pi-scaled"
trig functions allow for pathological matches like
"tanpi(tanpi(tanpi(tanpi(1/tanpi(tanpi(ln(x))))))) = 1/3" mainly
because of the way the derivative scales up with each nested call of a
trig function. Al address this, I limit the ranges of sin, cos and tan
to +-pi (after applying the scale factor). Of course you don't want to
reduce symbol weights so drastically, and more realistic weight
adjustments (like "ries 2.5063 --symbol-weights 8:T") are more useful.
 Due to the argument syntax you cannot set a negative weight, and
large weight values automatically disable the symbol.
 Currently, The automatic setting of k.vanished_dx for large values
causes tautology errors when RIES is given large values. For example,
currently "ries -l3 1234567890" gives "x^2-sqrt(x)^4 = 2^8"
20121210 Pass parameters using_x, a_minw, a_maxw, etc. through the
recursive expression generation routines, rather than having them be
global variables.
 Add tautology warning messages when the target is very large or the
--derivative.margin is too small. Here are examples, in all cases using
a command like "ries -l4 4243743 --derivative.margin 2.0e-11", RIES
gives an answer like "x^2-sqrt(x)^4 = 1/-(2^8)"
 Target min. --d-m
 42437432.1 10.0e-9
 4243743.1 188e-11
 424374.1 11.8e-11
 112474.1 3.0e-11
20121211 Add a bunch of comments in ge.2(), etc.
 Add cv.simplify rule [Arr]=[B] -> [A]=[B]
20121212 Fix a bug in infix.1(): "7 times negative x" is now displayed
as "7*-x", before this fix it was "7 -x" which looks too much like
"7-x"..
20121214 Improve k_vanished.dx checks and error messages.
 It has been a while since I did benchmarks using the standard test
(2.506+pi/10^4 with no options), so I am going to run the benchmark
on all my machines again.
 Although I have been doing the same test since the beginning (see
above under the dates 200003xx, 200310xx, 20050715, 20090216 and
20090515) the benchmarks are not directly comparable because of many
changes in the code. Improvements in error-checking, pruning rules,
the addition of the tangent function, and recent fixes to the
derivatives of trig functions have all affected the number of
generated expressions and therefore the running time and memory usage
for a given search level. (For details of all the changes, see above
under the dates 20101218, 20111219, 20111220, 20111221, 20111226,
20111228, 20111229, 20111230, 20120505, 20120613, and 20121208). Here
is a representative sampling of statistics and running times from
"ries 2.5063141592653589 -l4" over the past few years:
 Date and time expressions complexity time (Xeon
 YYYYMMDD.hh:mm LHS RHS LHS RHS equations memory E5520)
 20090217.02:36 1344958 1080143 76 70 1.453e12 155.3M 10.6s
 20090515.18:16 1344952 1080143 76 70 1.453e12 155.3M 10.5s
 20090516.01:13 1332922 1360355 76 71 1.813e12 172.4M 12.6s
 20111219.21:00 1332922 1360355 76 71 1.813e12 172.4M 12.2s
 20111220.23:28 1327806 1360355 76 71 1.806e12 172.1M 12.4s
 20111221.23:46 1327806 1360355 76 71 1.806e12 172.1M 12.3s
 20111226.18:27 1327757 1360355 76 71 1.806e12 172.1M 12.4s
 20111228.22:24 1251691 1280315 76 71 1.603e12 162.1M 11.60s
 20111230.06:26 1020295 1252891 76 72 1.278e12 145.5M 11.61s
 20120505.22:54 1032922 1281451 76 72 1.324e12 148.2M 11.55s
 20120613.23:35 1032922 1281451 76 72 1.324e12 148.2M 11.82s
 20121209.09:58 1032947 1281451 76 72 1.324e12 148.2M 11.80s
 20121210.07:24 930662 1127640 76 72 1.049e12 131.8M 10.74s
Here are the current times for each level from -l0 to -l6, on all
machines back to the 800 MHz G4 (which is 10 years old).
 PowerPC PowerPC Core 2 Xeon Core i7
 G4* G5* Duo E5520& 2720QM#
 Level -mem- eqns. 800 MHz 2.0GHz 2.23GHz 2.27Ghz 2.20GHz
 -l0 1.11Mb 6.89e7 0.333s 0.085s 0.035s 0.026s 0.022s
 -l1 3.47Mb 7.18e8 0.954s 0.304s 0.125s 0.102s 0.077s
 -l2 11.9Mb 8.46e9 3.495s 1.245s 0.518s 0.436s 0.364s
 -l3 39.3Mb 9.23e10 14.70s 6.407s 2.724s 2.125s 1.507s
 -l4 131 Mb 1.04e12 69.29s 34.13s 12.75s 10.70s 7.760s
 -l5 441 Mb 1.19e13 % 165.0s 64.15s 51.59s 38.29s
 -l6 1.51Gb 1.39e14 794.0s 322.2s 252.3s 190.1s
 * For the G4 and G5 I was using a 32-bit binary, so the memory usage
 was a bit smaller.
 % The G4 did not have enough memory to run the -l5 or -l6 tests.
 & The Xeon E5520 has a base speed of 2.27Ghz, but during these tests
 all cores were otherwise idle, so it was running near its single-core
 turbo speed of 2.53 GHz.
 # The i7-2720QM has a base speed of 2.2 GHz, and all cores were otherwise
 idle during these tests. However, its turbo frequency depends on the
 amount of recent activity. Each test was run after the CPU had been
 completely idle for 1 minute. For tests under 25 seconds, this means
 the clock speed was near the maximum 3.3 GHz.
20121215 Standalone math source (msal_math64.c) now includes SIN and
COS functions.
20121216 Fix a simple bug in symstrcmp that prevented
--find.expression from working.
20121217 Add --explicit-multiply option.
20121218 Allocate twice as much space for the RHS part of the equation in
try.solve and in report.match
20121223 Slight optimization inside the next symbol test loop of ge.2.
 Fix a few compilation warnings (missing prototypes; mismatched
integer types in some debug printfs)
20130121 Fix bug in -S option when used with --numeric.anagram (that
made it act as though the -S option had not been given).
20130127 While working on a simple Perl script to generate OEIS sequence
A005245 I discovered that the command:
 ries --one-sided -ie -S1+* 23 -l5
does not find any answer. I suspect this is a bug in complexity limit
optimizations. {After investigation I figure out it's simply that
MAX_ELEN is too low.}
20130129 Do a thorough inspection of all places where I check
expression lengths against MAX_ELEN and find a few places it can be
improved to make maximum use of the avaiable space. Allocate one more
symbol for certain expression scratch buffers. It should be able to
use all MAX_ELEN symbols now. Add RIES_MAX_EXPRESSION_LENGTH option, which
you can give when compiling RIES to make a RIES that handles longer
expressions.
 Make "exhaustion timeout" message suggestions more relevant;
don't show it at all if a result has been given; don't show "(for more
results, use the option '-l3')" suggestion if exhaustion timeout error
was given. Also add NO_IDENTITY_OPTIMIZATION code (currently disabled).
20130130 Add a 32-bit signed integer data type and corresponding
ifdefs and runtime sizeof testing. Symbol attributes mark is now 32
bits. Add rules for associative operators: for example, it now prunes
[ABC++] = A+(B+C) in favor of [AB+C+] = (A+B)+C.
20130201 Now using debug_B for the derivative prune messages formerly
displayed by debug_e and debug_E; re-use debug_E for "rejected
(duplicate value)"; debug_G now (properly) prints its message only
when an expression has actually been added, and uses infix notation
because -DG/-Dg may be useful to users who wish to use RIES to
generate a database of expressions for processing by another program.
20130203 Add Lambert W function define and test loop; then move the
test code to msal_math64.c. To see LambertW tests, compile with
RIES_USE_SA_M64 defined then invoke RIES with option -Dz.
20130218 "--max-match.distance 0" now causes RIES to exit if and when
it gets an 'exact' match. Add a roundoff-error disclaimer about
'exact' matches (unless -i was given).
20130228 Add null argument '-' (useful if you want to end a sequence
of arguments to something like --symbol.weights).
 -S, -O and -N now pre-empt each other in the order they are given
(which is necessary for users who combine different profiles).
 Add the -E option (replacing the old "or-mode" behavior of -S; this
will be necessary for future optional functions like A and W, and also
facilitates combining multiple profiles.
20130301 Change almost all occurrences of 'double' declarations to
one of
 ries_val: The value of an expression or subexpression
 ries_dif: The value of a derivative, error, uncertainty, the distance of
 a match, etc.
 stats_count: For counting generated expressions, equations, etc.
 time_flt: For measuring time and memory usage
 Add init.formats()
20130302 Move Gamma function code (now tested) to msal_math64.c. As
with LambertW (see 20130203) it is not yet available in expressions or
equations but you can see the internal tests of Gamma by compiing with
-DRIES_USE_SA_M64, then invoke RIES with the option -Dz.
 Add default --include/-p profile feature. It now looks in
getenv("HOME") or getenv("USERPROFILE") in Windows for either
"ries_profile.txt" or ".ries_profile" (either name works in all OS's).
This profile is loaded first before any of the other options; however
if you give no options RIES still prints the brief help and exits.
 Add a '-p' option that:
 - If given first, disables the default profile
 - If given anywhere else, loads the profile at that point (so you can
 give an option which is applied before loading the profile by
 e.g. "ries -p foo bar -p baz qux")
 Make the '--version' option display the path to the defaults file,
to aid field diagnostics.
20130303 Add attribute tags (TAG_INT, TAG_RAT, TAG_NONE) and implement
their calculation. This provides a more robust way to implement the -i
and -ie options, and more importantly it paves the way for new
features.
 Add --absolute-roots as synonym for -x; add --relative-roots option
(opposite of '-x').
 Use newly-added attribute tags to implement -i option (replacing
many calls to float() with hopefully less aggregate work; look in
exec() and search for 'g_restrict.subexpr'.)
20130305 Add --no-solve-for-x option (which merely undoes
--try-solve-for-x). Add setup_abc_mmw() as part of refactoring the
handling of restricted symbolsets.
20130306 Add --symbol.names option; this involves several changes to
how we define the symbols, all of which will help with future
improvements like user-defined constants.
20130307 '=' symbol can now be renamed. Add -r and -re options to
ensure that all solutions are rational when solved for x. (These are
just shorthand for using -N to exclude lots of symbols). Increase
allocation of temp buffers in infix.1, fixing a crash that happened
when using an all-seft-b symbolset and solving for x, e.g. "ries
3.1415926535897932 -SeElq -s".
20130308 Add --no-refinement option, which causes it to never
decrease best.match, and therefore print *all* matches that do better
then the specified distance. This produces results similar to the ISC
(except still ordered roughly by increasing complexity) and may be
useful to folks (like the 137 cultists) who are happy with any formula
within a known error bound. It required adding sym.strneq,
symstrsymstr, unique.eqn and the g.matches memory block (none of which
are used unless you choose the option).
 Add the --max-matches option, 100 by default.
20130309 Add -n as a synonym for --max-matches, --integer-subexpressions
as a synonym for -i, and --rational-subexpressions as a synonym for -r.
20130310 Add the g_nr_deltas array to provide another type of pruning
for equivalent answers for --no-refinement.
 Add the "-c" (alias --constructible-subexpressions) option which is
similar to -r but also allows phi, square and square root. Do a little
refactoring of how the -D, -E, -i, -N, -O and -r options are handled.
 Begin refactoring the tags manipulation to use TYPE_xxx values: since
each attribute is a subset of another (e.g. the integers are a subset
of the rationals) I don't really need to use bit-fields.
20130311 Continued refactoring of tags, add TYPE_CONS and TYPE_ALG
tagging. Implement -a option; add tagname().
20130312 Add --rational-exponents option, which tests the tag of the
argument and generates the new error ERR_EXEC_ILLEGAL_EXPONENT; the
error is also generated if the argument contains x on the assumption
that the user probably does not believe x to be rational if he is
looking for an answer in the form of an algebraic number.
 Add --rational-trig-args option, which works similarly and also sets
k_sincos_arg_scale to pi; it generates the new error
ERR_EXEC_TRIG_ARGTYPE.
 Add -a (algebraic subexpressions) option, using these new restrict
options. To support this, the target number is automatically tagged
with the same type as the selected restriction (-r, -c, -a) unless it
is obviously an integer or half- or quarter-integer; if no restriction
option was given we tag x as transcendental. There are now four
classes of restrictions with easy command-line options: -i, -r, -c,
and -a.
 Add ":.:." syntax for defining a blank-space character in symbol
names.
20130314 Add --any-exponents and --any-trig-args options.
20130317 Compute rv_maxint and disallow -i when target is too large;
remove utf8-related code (obviated by --symbol.names)
20130318 Expand brief.help.
 Add "-l" option (without a numeric argument, i.e. a bare "-l" rather
than something like "-l3") to restrict answers to those that have
Liouvillian numbers as roots. With -l it finds that 1.632526919438153
is sqrt(2)^sqrt(2) and finds that 1.132997565885066 is a root of
x^5+x=3, but does not find that 2.31645495878561 is the root of x^x=7.
Along with the full unrestricted defaults, and the option "-Ox", this
makes for a rather full set of options to select popular classes of
numbers as the roots of the reported results.
20130219 Test g_restrict_exponents in functions 'l', 'E' and 'L'.
20130613 Add --no-slow-messages option; auto-set k_min_best.match and
k_vanished_dx for small targets.
20130626 Use "%Lg" in various debug printfs to support RIES_VAL_LDBL.
Add macros EXP, FABS, etc. and (ries_dif) typecasts in several places;
get long double precision pretty much working.
20130801 Initialize k_ulp sooner so it can be used by -ce
 --max-match-distance and --match-all-digits options now cancel each
other.
20130803 In --symbol-names, allow redefining a symbol to itself (e.g.
':^:^' when -F format is selected) by accepting space_sym definition
only once.
20130805 Benchmarks of standalone maths library vs. standard libm.
Using the command "time ./ries-libm 2.5063141592653589 -l5
--max-match-distance 1e-10", and SIMULTANEOUSLY running the same
command with the "ries-sa-math" library, the time is 0m45.373s with
the sa-math library and 0m45.568s using libm (the sa-math library is
actually faster).
20130809 Add ieee.paranoia()
20130810 Alphabetize the order of sections in parse_args()
20130811 Add ERR_EXEC_ILLEGAL_DERIV; exec() checks for overflow and
NaN in derivative calculations in several operators
 Print error and exit if target value is zero.
20130812 Pass root directly to report.match to get full precision in
the case where ries_val has much more precision than ries_dif.
Increase precision of calculations of pi and e in init.numerics. "long
double" now gives 31 usable digits on the PowerPC G5 (where GCC 4.0.1
and later provide double-double arithmetic), and show.version now
shows the architecture (PPC/Intel) and precision:
 ries --version
 ries version of 2013 Aug 13, Copyright (C) 2000-2013 Robert P. Munafo
 architecture: PowerPC
 precision: long double (33 nominal, 31 usable)
 mathlib: standard
 profile: -p/Users/munafo/.ries_profile
 RIES is free software; see the source for copying conditions. There is NO
 warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 Dual-G5 /Users/munafo
 : ries 2.50618 -x
 Your target value: T = 2.50618 mrob.com/ries
 x = 5/2 for x = 2.5 {50}
 x = e^3/8 for x = 2.510692115398458467616066206823 {67}
 x = sqrt(2*pi) for x = 2.506628274631000502415765284811 {55}
 x = x"/(1+9) for x = 2.506184145588769256292940922378 {70}
 x = pi^(ln(sqrt(6))^2) for x = 2.506182746702861729338863939787 {84}
 x = e^(5"/sqrt(3/7)) for x = 2.506182083318281693444031686493 {99}
 x = ln(pi*e+1+e) for x = 2.506180094504077068176274103417 {96}
 x = e^(log_(pi^2)(7))+1/6 for x = 2.506180001393763345666526465285 {107}
 x = sqrt(1/((1/e-ln(6))+sqrt(x)))
 for x = 2.506180000397503542091976573633 {118}
 (for more results, use the option '-l3')
 log_A(B) = logarithm to base A of B = ln(B) / ln(A) sqrt(x) = square root
 e = base of natural logarithms, 2.71828... A"/B = Ath root of B
 ln(x) = natural logarithm or log base e pi = 3.14159...
 --LHS-- --RHS-- -Total-
 max complexity: 67 62 129
 dead-ends: 2842183 5668638 8510821 Time: 2.106
 expressions: 198545 370110 568655
 distinct: 103576 107460 211036 Memory: 13248KiB
 Total equations tested: 11130276960 (1.113e+10)
By comparison, typical run time for the same command on PMG5 in double
precision is about 1.49, and memory usage is 10624KiB.
20130816 Add ries_strncpy, ries_intpow, ries_to_digits,
ries_snprinf_int, ries_strlen, ries_spfg, spfg and spff. Use these
in a few debug statments.
20130818 Remove trailing 0 digits from ries_spfg output (for example:
2.50618e+12 instead of 2.50618000e+12); fix a few compilation
warnings.
20130820 Add ries_spfg_test and msal_test_spfg
20140226 Eliminate conversion to double in ries_to_digits, fixing an
'exponent adjust failed' error when trying to print values like
1.23e-789 (with long double which supports exponents down to e-4951)
20140227 Minor refactoring; add did_newton parameter to report_match
20140228 Add some logic to check.exact_match, duplicating some of the tests
done by check.match.
20140303 Generate ERR_EXEC_ZERO_DERIV error in '^', 'L' and 'v'
operations when the derivative calculation underflows to zero. This
eliminates several recently-reported tautology errors.
20140831 Expand brief_help() a little,
20141014 Add --min-equate-value and --max-equate-value options.
 20141105 Add Lambert W function (duplicating and expanding on work by
Mark Shoulson).
 20141106 Increase weight of 'tan' (that is, the tangent function when
--trig-argument-scale is not the default) to 6, bringing it in line with
'tanpi' (where the weight was always 6).
 Bug-fixes in Lambert derivative and error handling, in debug_z mode
call new test routine msal_test_lambertl
 20141117 Make -i (integer subexpressions) more efficient by restricting
symbolset as is done for -r. All class restrictions now also exclude W.
 20141122 RIES license is now GPL version 3 (formerly GPL v2). Reorder
the info in --version a little bit.
 20141207 Add test of each subexpression against g_restrict.subexpr so
that if you do e.g. "-a -EL" the L will be enabled but effectively
ignored; however I also canonize the unofficially supported "-a -Ep"
by mutating the types of the predefined constants in a similar way as
was already done for the target value. This is how user-defined
constants without a user-defined class will need to be handled.
 Add --any-subexpressions to enable the old behaviour, e.g. the user
can use "-a --any-subexpressions -EL" to get the equivalent of the old
"-a -EL".
 20141212 struct form now includes stk[] and arg1[], supporting rules
that check the first argument of binary operators; add one rule using
AM_a1_e.
 20141213 Add AM_a1_1 rules; #qualify# tests show that several more
results are found in -l5 tests, indicating that these rules increase
pruning for deep searches.
 20141216 Replace all "sprintf" with snprintf.
 20141217 Add #defines for most of the snprintf tempbuf sizes.
20160104 Enable cube root of a negative argument ("-9.8222414378011"
vector in qualify tests it)
 20160131 Add mem_used_bytes for more precise reporting of how much
memory would be needed for small tasks; to support this I did some
signed-vs.-unsigned cleanup.
 20160423 Finish re-indenting the big block of options tests in 
parse.args()
20170211 Add '--show-work' as a synonym for '-Ds'
20180713 Add '--max-trig-argument' option.
 20180802 Rename '--max-trig-argument' to '--max-trig-cycles'
20210417 Clean up arctan code and get it working.
 20210418 Better formatting of atan2 in infix.1(), don't print
2nd argument if it is 1.
20240422 Legend now explains the two forms of atan2 (one-argument
and two-argument) separately, giving only one explanation when only
one of the forms appears in output.
BUGS and TO-DO
 See "UNFINISHED WORK" section above
THE CHECK_SIDES PARADOX
The RIES algorithm maintains a single list of RHS and LHS expressions
sorted in numerical order. These are distributed pretty much randomly,
and when a new expression is added, RIES checks the preceding and following
list items to see if it can form an equation with the new item.
 The closeness of a match is the difference in values divided by the
derivative of the LHS:
 closeness = |LHS(x)-RHS|/(d/dx LHS(x))
If the new node is an LHS, then the values of LHS(x) and d/dx LHS(x)
are the same for every match-comparison that is made. One RHS in each
direction (upward and downward) is all that needs to be tested,
because any further RHS's will generate a larger value of
|LHS(x)-RHS|.
 However, when the new node is an RHS, the value of d/dx LHS(x) will
differ for each LHS that is found while scanning upward and downward
for candidate matches. Therefore, even after finding a new optimal
solution, the possibility exists that there might be another
even-closer optimal solution if you keep scanning further, coming from
an LHS that has a much lower value of d/dx LHS(x).
 Here is a concrete example. La and Lb are two LHS's with derivatives
of 1 and 10. Ra and Rb are two RHS's. They are shown here as if laid
out on a number line to make the example clearer. The nodes are
inserted in the order: La, Lb, Ra, Rb.
 value: 1 2 3 4 5 6 7 8 9
 Ra Rb La Lb
 d/dx: - - 1 10
When Ra is inserted, the match La=Ra is found, with closeness
(8-1)/1=7. Then Rb is inserted, and a new match La=Rb is found, with
closeness (8-7)/1=1. These are the only two matches it will report.
 However, the match Lb=Ra is has a closeness of (9-1)/10=0.8, which is
closer, and Lb=Rb has a closeness of (9-7)/10=0.2 which is closer
still. Lb=Ra should have been reported instead of La=Rb. The program
should report La=Ra, Lb=Ra, and Lb=Rb (in that order).
For a real-life example, use the command:
 ries .328106566874978253 -l-4 --trig-argument-scale 1 -NT -Dy0
which only gives one result, "5 x = sqrt(e)". With -D0 it dumps the
entire table of values on each complexity pass. Look through this
output for the first occurrance of [xrS] and [9r] together:
 28 xrS { 35} = 0.093664890868384448 , dx = 9.2481888991224963
 29 xp/ { 34} = 0.10443956395812863 , dx = 0.31830988618379069
 30 xs { 24} = 0.10765391922648457 , dx = 0.65621313374995649
 31 x3/ { 35} = 0.10936885562499275 , dx = 0.33333333333333331
 32 9r { 26} = 0.1111111111111111
 33 8r { 26} = 0.125
The output shows that [9r] is inserted a few passes later than 'xrS',
so check.sides is looking at [9r]'s neighbors [x3/] and [8r].
 For a while in 201112 I made RIES look past the nearest neighbor,
and this command reported the result "sin(1/x) = 1/9". But as you can
see, sin(1/x) matches 1/9 more closely than x/pi, x^2, and x/3,
because the derivative of sin(1/x) is so much higher.
This is an example of the above-described problem, and is the reason
why for a while RIES was checking multiple neighbors on each insert.
The Paradox:
Despite the foregoing, the present algorithm (in which check.sides
only looks at the one closest neighbor on both sides of a newly-added
expression) turns out to work very well.
To understand why, look at the output of ries 1.506591651 -Dy0 and
find the first appearance of [1p6*-] :
 448 x1+p^n { 56} = -17.937341258252982 , dx = -22.481451854903934
 449 1p6*- { 51} = -17.849555921538759
 450 xe+sn { 51} = -17.84955591743638 , dx = -8.4497469589180909
 451 xTe/ { 49} = -17.762306262640966 , dx = 857.98450178681276
 452 x7^n { 46} = -17.61849853924631 , dx = -81.859931782354067
When [1p6*-] is inserted, the LHS expressions [xe+sn] and [xTe/] are
already present. [xTe/] has a much bigger derivative (over 100 times
as large), so it looks like a good candidate for a match that would be
missed if check.sides only looked at the first neighbor of [1p6*-].
However, the distance in x values from [1p6*-] to [xTe/] is over 20
million times larger in magnitude than the distance from [1p6*-] to
[xe+sn], so the higher derivative of [eTe/] doesn't stand a chance.
When checking for other similar cases, the same thing always happens:
when an RHS and LHS forms a new record close match, any other LHS's
in the area don't come anywhere close to being another new match.
The reason for this is in the statistics. Going back to the
.328106566874978253 example, consider the range of values between 0
and 1. By the time there are 1000 expressions in this range, the
average distance between expressions will be about 0.001. However, the
distance between the two *closest* expressions will be much smaller,
somewhere on the order of e/10^-6.
 Now consider what happens when you insert another 1000 expressions
at random places in the range (0..1). There is a reasonably good
chance that one of these new 1000 points will come closer to an
existing point than any of the other old points was. This will be a
new match. However, the odds of having *another* new match at the same
time are very very low -- about 1 in 1000. In other words, in order
to get a situation where there are two good LHS matches near an RHS,
all three have to be within 10^-6 of each other.