/* ries.c RIES -- Find Algebraic Equations, Given Their Solution Copyright (C) 2000-2023 Robert P. Munafo This is the 2023 Aug 01 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:
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 ...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