konn / computational-algebra

General-Purpose Computer Algebra System as an EDSL in Haskell
http://konn.github.io/computational-algebra/
BSD 3-Clause "New" or "Revised" License
92 stars 9 forks source link

Floating-point algorithms? #14

Open sheaf opened 3 years ago

sheaf commented 3 years ago

Apologies if I've missed something, but I couldn't figure out how to use solveM or solve' with floating point coefficients (e.g. Double). Is this possible?

In my situation, the algorithms using Fraction Integer seem much too slow. Something that takes less than a second in Mathematica or Macaulay2 doesn't even terminate after 15 minutes with solveM/solve'.

Here's an example of a system that I've been working with (I cleared the denominators manually, for simplicity)

p = -495 - 5625 * t + 27891 * t^2 - 33372 * t^3 + 13824 * t^4 - 2367 * t^5 + 6345 * u + 25110 * t * u - 170721 * t^2 * u + 230166 * t^3 * u - 110682 * t^4 * u + 15354 * t^5 * u - 18765 * u^2 + 22140 * t * u^2 + 148581 * t^2 * u^2 - 256878 * t^3 * u^2 + 148446 * t^4 * u^2 - 24300 * t^5 * u^2 + 25110 * u^3 - 100440 * t * u^3 + 56970 * t^2 * u^3 + 4860 * t^3 * u^3 - 29160 * t^4 * u^3 + 16740 * t^5 * u^3 - 12960 * u^4 + 51840 * t * u^4 - 27540 * t^2 * u^4 + 33750 * t^3 * u^4 - 21600 * t^4 * u^4 - 8370 * t^5 * u^4 - 3240 * t^2 * u^5 - 25920 * t^3 * u^5 + 25920 * t^4 * u^5

q1 = 52785 + 56943 * t - 4037688 * t^2 + 15883803 * t^3 - 26239221 * t^4 + 20703357 * t^5 - 7724700 * t^6 + 989361 * t^7 - 66015 * u + 487134 * t * u + 27385938 * t^2 * u - 113650722 * t^3 * u + 181819755 * t^4 * u - 135044982 * t^5 * u + 48046662 * t^6 * u - 6530274 * t^7 * u - 2478195 * u^2 + 6294834 * t * u^2 - 111627882 * t^2 * u^2 + 417376962 * t^3 * u^2 - 654613569 * t^4 * u^2 + 478520946 * t^5 * u^2 - 162633582 * t^6 * u^2 + 21861414 * t^7 * u^2 + 13768380 * u^3 - 50360616 * t * u^3 + 247545072 * t^2 * u^3 - 696128256 * t^3 * u^3 + 1052548344 * t^4 * u^3 - 789620400 * t^5 * u^3 + 272718252 * t^6 * u^3 - 36077400 * t^7 * u^3 - 28987470 * u^4 + 118189368 * t * u^4 - 321295572 * t^2 * u^4 + 526662486 * t^3 * u^4 - 679607010 * t^4 * u^4 + 541760724 * t^5 * u^4 - 204607620 * t^6 * u^4 + 26871750 * t^7 * u^4 + 27060480 * u^5 - 113821200 * t * u^5 + 244827360 * t^2 * u^5 - 186998220 * t^3 * u^5 + 124853400 * t^4 * u^5 - 150553080 * t^5 * u^5 + 83446200 * t^6 * u^5 - 9491580 * t^7 * u^5 - 9331200 * u^6 + 39074400 * t * u^6 - 91154160 * t^2 * u^6 + 40979520 * t^3 * u^6 - 5093280 * t^4 * u^6 + 57678480 * t^5 * u^6 - 44634240 * t^6 * u^6 + 3013200 * t^7 * u^6 + 233280 * t * u^7 + 3265920 * t^2 * u^7 + 14929920 * t^3 * u^7 - 22161600 * t^4 * u^7 - 3732480 * t^5 * u^7 + 9331200 * t^6 * u^7

q2 = 919350 - 1206090 * t - 18350955 * t^2 + 58802085 * t^3 - 67806099 * t^4 + 35782614 * t^5 - 8371728 * t^6 + 651591 * t^7 - 5576850 * u + 32915970 * t * u - 20491785 * t^2 * u - 89979930 * t^3 * u + 154823643 * t^4 * u - 98196624 * t^5 * u + 26236818 * t^6 * u - 2518830 * t^7 * u + 9173250 * u^2 - 91559970 * t * u^2 + 183517245 * t^2 * u^2 - 118728180 * t^3 * u^2 - 7137315 * t^4 * u^2 + 55358640 * t^5 * u^2 - 33001290 * t^6 * u^2 + 5753700 * t^7 * u^2 + 1069200 * u^3 + 54276480 * t * u^3 - 129054870 * t^2 * u^3 + 117919800 * t^3 * u^3 - 53146530 * t^4 * u^3 - 19948680 * t^5 * u^3 + 38889720 * t^6 * u^3 - 9797760 * t^7 * u^3 - 11372400 * u^4 + 12733200 * t * u^4 - 9379800 * t^2 * u^4 - 12830400 * t^3 * u^4 + 49223700 * t^4 * u^4 - 15870330 * t^5 * u^4 - 16692480 * t^6 * u^4 + 7690950 * t^7 * u^4 + 4665600 * u^5 - 4665600 * t * u^5 + 21724200 * t^2 * u^5 - 38296800 * t^3 * u^5 + 5297400 * t^4 * u^5 + 1628100 * t^5 * u^5 - 356400 * t^6 * u^5 - 2762100 * t^7 * u^5 - 583200 * t * u^6 - 11664000 * t^2 * u^6 + 23328000 * t^3 * u^6 - 4422600 * t^4 * u^6 + 6536700 * t^5 * u^6 - 2721600 * t^6 * u^6 + 753300 * t^7 * u^6 - 874800 * t^4 * u^7 - 4665600 * t^5 * u^7 + 2332800 * t^6 * u^7

These are polynomials in two variables t, u. Both Mathematica and Macaulay2 can solve the system {p, q1, q2} in less than a second, whereas neither solveM nor solve' can finish under 15 minutes.

In this case, I'm looking for the real solutions with t,u in the unit interval, which are:

polynomial_system

I also tried to use Fraction Int, but a lack of instances prevented that from working.