Interpolation over an integral domain needs proof

[wbhart]

Currently we interpolate polynomials over an integral domain by first converting to Newton basis, doing an interpolation by the method of divided differences, then converting back to monomial basis. The problem is, if there is no polynomial over the integral domain which interpolates the initial values, the interpolation stage will do an inexact division. Currently we raise an exception if this happens. But what we haven't proved is that if there is a polynomial over the integral domain that interpolates the input values, then no fractions can occur at this point. This isn't trivial to prove (usually textbooks assume all the differences are units to make the proof work). However, it does seem to be true, both experimentally, and using a symbolic ring as input. What is missing is a mathematical proof that it always works without going to the fraction field internally. Beware that there is no obvious proof, and one must check carefully that any purported proof really checks that no inexact division could occur in the case of an exact interpolation, and that the same proof doesn't show that inexact divisions never happen, even when there is no polynomial interpolating the inputs!! The difficulty in the proof is that fractions could cancel out in the arithmetic in the inner loop of the interpolation. Working backwards naively from a valid interpolated polynomial isn't enough to show that this didn't happen.

[wbhart]

I am now occasionally hitting cases where the division is not exact, even when the interpolation can be done. I think it should be easy enough to generate examples since I've seen it randomly twice now. It's always Flint that says the division was inexact.

[wbhart]

Here is the counterexample:

x = Nemo.fmpz[-19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] y = Nemo.fmpz[105376795671912660390950794733193591478185572967687766712482771047560111844653631916205711049033, 12246053179864692206244836485475126466736376418650194736205466830505588580059416275520114862373, 1259203390208420569633553017292286248955410640987637614780340054562893412864553709252773553195, 112883232062918232174923890175187220899286614152347953291828052584483704976289537426697906181, 8668482199227025131256531735717138211259041662134612570355826725929935396496775078402374005, 558214884177553770313195091637122291610293950165013955245298602100864119114980000513046053, 29370851401402266090340124017345188788533230709235022861310786268501317991425074567441063, 1222602938121337891538787819891462621011852001151367397353942174061793295214225122592965, 38657409685493755237258333010628917423647082033962623342536627318819370491099330531521, 881011704478313692815185634396762152272070231103089652312881734198252670350958254245, 13503107090107449958413773918490034035811331960708431425456256884604866310225401923, 126630608432512492343026933876936383646036801265905948731450325013969997606798853, 635239369893238466716907362805088515456023444033538925776177806924641176490285, 1391994300963597134040355543629067293743584983790701956381003812736577053861, 951789302196091257982280311846098129914762490522348662191021935121339135, 100875759778261026040583968796129264628484389091128245718162640381893, -729099670657846511246047687666220514494312822536384758504072007, -769603331839146317492859383801412708076340043178461829595, 1458970778809069326337537299376030882639343327451, -7167982347705640439764333055704018204945925, 299372439593338888601285820783474092205269626213, -504592550638244504420806507215373988920768434704708034267, 455480913586038299368470363235361809592458746289176697702383575, 58747384244885713513659364491937888815976256034328111346847604281541, 451882843334833437916996161700513728895825120080678007879268128617026865, 672242968035942584879176959920558695128646391832791482272038487163942900133, 323679479015376548549412312973014442608198860712577522558618296671092351694323, 68285510557492315744063773385257933763083507498471285399958999037541737962455045, 7670026961397366932731033198187653449899212098281904420730851648783322048208420381, 523912168204705785980378726200423239383152182412863542889267318996372877115636919205, 23927136414631216976494012790928582672954730205223167444635151814718761626184925593903, 783619543556484278851402912493070814108663887059497117566405863833902234743242348291013, 19409180905260096762034048185154361596736812978161081724782197297327860923599150670785065, 378935780365528590330748733808877645656880714087513261192600222101197848975435157488573221, 6026043798215065661329042194369612501371341215607917757503892720984950105512643895080845195, 80149652296292166883452746135295215910706530955303060716383115947795121028319378044717472773, 911127664246654392133672616805069950489783840408911037977421161556050633621530950770411828053, 9012785486630311648044488891610441365279196037665865412371688528950425725560293641628715598885, 78753839906119380429576884732296055019298223957112598076798589738901626643333295411688691394311, 615663959921317343463947166298756615497885690746797932681508956024918160644419916815494218662085, 4353042757813749034130860179719753495875056468252525209197643063088277702855443593819553945176993]

[wbhart]

Actually, this isn't a counterexample after all. If I interpolate over Q, the resulting poly doesn't have coefficients over Z. This came from a failing case in our solve_poly benchmark, which indicates that there is a bug somewhere else.

[thofma]

Do you have the failing input of solve_poly?

[wbhart]

Sorry, I forgot to update the ticket. The reason for the failure is that Flint only computes the determinant up to the sign. I have fixed this in Flint by explicitly computing the sign of the determinant. We will need to do the same thing in Nemo, as soon as I can figure out what the rabbit warren of solve_rational and solve_with_det functions is actually supposed to do.

[wbhart]

We may be able to prove this by staring at this [1], as interpolation is equivalent to solving

(vandermonde)*(coeffs) = (values)

The method of divided differences (which we are using) looks like it boils down to this.

https://proofwiki.org/wiki/Inverse_of_Vandermonde%27s_Matrix

[wbhart]

The proof is given for integral domains on page 383 of the following paper [1]: "Chinese remainder and interpolation algorithms", J. D. Lipson 1971.

[1] https://www.google.com/url?q=https%3A%2F%2Fdl.acm.org%2Fcitation.cfm%3Fdoid%3D800204.806309&sa=D&sntz=1&usg=AFQjCNFy55xr2eiiysANAf1BDYa7Z7oEzg