-
Notifications
You must be signed in to change notification settings - Fork 0
/
ries.c
13850 lines (12380 loc) · 490 KB
/
ries.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* ries.c
RIES -- Find Algebraic Equations, Given Their Solution
Copyright (C) 2000-2022 Robert P. Munafo
This is the 2022 Jun 30 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 -o ries ries.c -lm
If the compilation fails and reports errors like "undefined reference",
try moving the pieces around: "gcc ries.c -lm -o ries". 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_GSL
Use this option to link to GSL, the GNU Scientific Library, which
provides many more special transcendental functions. Information on
the library can be found at https://www.gnu.org/software/gsl/
-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
@ (a b -- b a) Swap top two stack elements
*/ /*
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