Optimal solution to find the coefficients from an unknown polynomial function

Question

Suppose I have a polynomial f(x) where all the coefficients a(i) are integers in the form of a function f, implement a function find_polynomial that takes f and returns the coefficients a(i) as a list of the form:

[a(0),a[1]....a[n]] (n<=15)

where a[n] is the last non-zero coefficient of f(x). The only other information about f(x) is that it is of finite degree (<=15). And the solution shall not import any module.

My approach was to try using the Lagrange polynomial, which is a more mathematical method of solving the solution, but it will not work for coefficient bigger than 10 ** 10), so I would like to know a much more stable solution which can solve all the unknown coefficients. Would like to see the passion for community to solve hard problems. Any help is welcomed :)

Unknown Functions:

def f(x):
    return 5 * x**3
def g(x):
    return 6 * x**6 - 2 * x**5 + 3*x + 5 * x**3
def h(x):
    return 5 + 3 * x**2 - 2 * x**11 + x + 3 * x**15
def z(x):
    return 5 + 3 * x**2 - 2 * x**11 + x + 3 * x**15-2*x**8
def p(x):
    return 5 + 3 * x**2 -31 * x**10 + x + 3 * x**15+13*x**8+14*x**9
def wrong(x):
    return 5 + 3 * x**2 -1 * x**10 + x + 10000000000 * x**15+13*x**8+14*x**9

Test Case:

print(find_polynomial(f))
print(find_polynomial(g))
print(find_polynomial(h))
print(find_polynomial(z))
print(find_polynomial(p))
print(find_polynomial(wrong))

Result:

(0, 0, 0, 5)
(0, 3, 0, 5, 0, -2, 6)
(5, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 3)
(5, 1, 3, 0, 0, 0, 0, 0, -2, 0, 0, -2, 0, 0, 0, 3)
(5, 1, 3, 0, 0, 0, 0, 0, 13, 14, -31, 0, 0, 0, 0, 3)
'''This ouput is wrong''' (5, 348713164801, -1133862589437, 1568627191296, -1243957041600, 638886422720, -226653467040, 57637291712, -10725815087, 1473646474, -149249101, 10998988, -573300, 20020, -420, 10000000004)

My approach:

def fact(num):
    if num==0:
        return 1
    else:
        return num*fact(num-1)
def poly(x):
    result=[]
    for i in range(x+1):
        y=fact(i)*fact(x-i)*((-1)**(x-i))
        result.append(y)
    return result
def find_polynomial(func):
    result=[0]*16
    def helper(len1,func,result):
        if len1==0 or func(1)==0:
            result[15-len1]=func(0)
            return result
        else:
            sum1=0
            polyx=poly(len1)
            for i in range(len1+1):
                j=(func(i)/polyx[i])
                sum1+=j
            result[15-len1]=int(round(sum1))
            return helper(len1-1,lambda x:func(x)-result[15-len1]*(x**len1),result)
    lst1=helper(15,func,result)
    def finalize(lst):
        if lst[0]!=0:
            return tuple(reversed(lst))
        else:
            return finalize(lst[1:])
    return finalize(lst1)

Do you have anything in particular in your code that you don't like (being not "python" enough), or are you asking for a general code/algorithm rewrite? — anatolyg, Jul 11 '18 at 09:03
I think this question might be better suited to the CodeReview stack exchange site https://codereview.stackexchange.com/. Usually working code that you want to improve is better received over there. — jambrothers, Jul 11 '18 at 09:30
I suspect that your mistake was to mark algorithmic problem with language tag. — MBo, Jul 11 '18 at 09:32
I’m unsure, reading your question, if the code work as intended. If it is, as @jambrothers suggests, it may be a better fit at Code Review. If not, please explain what differs between your actual results and the expected ones. — 301_Moved_Permanently, Jul 11 '18 at 09:40
@Mathias Ettinger@jambrothers, it is my bad to not point out, this solution will not work for coefficient larger than 10**10. Added that in the question. — Marcus.Aurelianus, Jul 11 '18 at 09:47
So you have numerical precision issues? There is approach: get `P(1)=Q`, then get `P(Q+1)=R` and express `R` value in base `Q+1`, but I'm afraid it suffers from the same problems for large coefficients. Is long arithmetics available? (I think Python has in-built multiprecision possibilities) — MBo, Jul 11 '18 at 09:58
@MBo, you can try it out in python or any language. I have done a similar solution before, it is more unstable than Lagrange polynomial. Python 3+ seems do not have long arithmetic, I think they got in-built multi precision. [Link here](https://stackoverflow.com/questions/538551/handling-very-large-numbers-in-python) — Marcus.Aurelianus, Jul 11 '18 at 10:00
Would appreciate if the community can support any thoughts....to solve the precision error in my code or provide an optimal solution for this problem. — Marcus.Aurelianus, Jul 11 '18 at 12:04
BTW, I missed that approach with `P(Q+1)=R` does not work for polynomials with negative coefficients — MBo, Jul 11 '18 at 17:42

MBo · Accepted Answer · 2018-07-12T04:01:10.863

I used algorithm POLCOE from chapter 3.5 of the book "Numerical Recipes in C", translated it to Python, then modified it for pure integer arithmetics to avoid rounding issues.

def f(x):
    return 1 - 2 * x + 3 * x**2 - 4 * x**3
def g(x):
    return 6 * x**6 - 2 * x**5 + 3*x + 5 * x**3
def h(x):
    return 5 + 3 * x**2 - 2 * x**11 + x + 3 * x**15
def z(x):
    return 5 + 3 * x**2 - 2 * x**11 + x + 3 * x**15-2*x**8
def p(x):
    return 5 + 3 * x**2 -31 * x**10 + x + 3 * x**15+13*x**8+14*x**9
def wrong(x):
    return 5 - 3 * x**2 -1 * x**10 + x + 10000000000 * x**15+13*x**8+14*x**9

def gcd(a,b):
    while b > 0:
        a, b = b, a % b
    return a

def lcm(a, b):
    return a * b // gcd(a, b)

def find_polynomial(func):
    n = 15
    m = (n + 1) // 2
    result = [0]*(n+1)
    s = [0]*(n+1)
    phi = [0]*(n+1)

    s[n] = m
    for i in range(1, n + 1):
         for j in range(n-i, n):
             s[j] -= (i - m) * s[j+1]
         s[n] -= (i - m)

    lc = 1
    for j in range(n+1):
         phi[j] = n + 1
         for k in range(n, 0, -1):
             phi[j] = k * s[k]+(j-m)*phi[j]
         lc = lcm(lc, phi[j])

    for j in range(n+1):
         ff = func(j-m)
         mul = lc // phi[j]
         b = 1
         for k in range(n, -1, -1):
             result[k] += b * ff * mul
             b = s[k] + (j-m) * b;

    for j in range(n+1):
         result[j] = result[j] // lc

    return result

print(find_polynomial(f))
print(find_polynomial(g))
print(find_polynomial(h))
print(find_polynomial(z))
print(find_polynomial(p))
print(find_polynomial(wrong))

result (note that I changed some your polynomials during debugging)

[1, -2, 3, -4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 3, 0, 5, 0, -2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[5, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 3]
[5, 1, 3, 0, 0, 0, 0, 0, -2, 0, 0, -2, 0, 0, 0, 3]
[5, 1, 3, 0, 0, 0, 0, 0, 13, 14, -31, 0, 0, 0, 0, 3]
[5, 1, -3, 0, 0, 0, 0, 0, 13, 14, -1, 0, 0, 0, 0, 10000000000]

for large power:

 def wronga(x):
    return 5 - 3 * x ** 2 -77777 * x ** 100 + x + 333333333 * x ** 15+13*x ** 8+14*x ** 9

 dont't forget to set proper max power:
 n = 100

 result:
 [5, 1, -3, 0, 0, 0, 0, 0, 13, 14, 0, 0, 0, 0, 0, 333333333, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, -77777]

Thanks! This approach is great for one extremely big coefficient. But if you try 2 relatively big coefficients. Like: def wrong(x): return 5 - 3 * x * * 2 -1 * x * * **100** + x + **100** * x * * 15+13*x * * 8+14*x * * 9 What I mean relatively big is these 2 100s. You can try it out in [here](https://repl.it/repls) for python IDE. — Marcus.Aurelianus, Jul 12 '18 at 00:52
The result seems not right for 2 relatively big coefficients. Like: def wrong(x): return 5 - 3 * x * * 2 -1 * x * * 100 + x + 100 * x * * 15+13*x * * 8+14*x * * 9, it is too long to post here. — Marcus.Aurelianus, Jul 12 '18 at 00:58
There is no big coefficients but big power. Perhaps you forgot to increase n to 100. It works well — MBo, Jul 12 '18 at 03:59
Seems it it worth to make function argument n/maxpower to provide minimal but sufficient work place. — MBo, Jul 12 '18 at 04:20
@ MBo, True, I did not notice that, trying to understand the formula you provide now. — Marcus.Aurelianus, Jul 12 '18 at 05:51

Optimal solution to find the coefficients from an unknown polynomial function

1 Answers1