0

I need to find Dp and Dq from this formula:

enter image description here

I tried to use ContFrac lib for this task:

pip install ContFrac

import contfrac

find_number = 2140e225
value = (1,math.sqrt(find_number))
conv = list(contfrac.convergents(value))
print(conv)

But the output is odd because number is too big:

0 + 1/(my_sqrt)

Method from SymPy also doesn't work and returns 0.

print(continued_fraction((1//sqrN)))

Article link

Update 1:

I was playing with SymPy and found that I should use this library sqrt implementation.

from sympy import sqrt

lm = continued_fraction((1 / sqrt(number)))
lz = flatten(lm)
print(list(continued_fraction_convergents(lz)))

But these will cause internal exception:

PrecisionExhausted: Try simplifying the input, using chop=True, or providing a higher maxn for evalf

But as I see I can't raise man for continued_fraction operation. And I'm also unable to pass number as decimal:

decimal.getcontext().prec = 125
lm = continued_fraction((1 / sqrt(decimal.Decimal(find_number))))

ValueError: expecting a rational or quadratic irrational, not 2.16152455423408e-125
Alexandr Dorofeev
  • 470
  • 1
  • 4
  • 12
  • I don't believe `ContFrac` supports this. Looking at the documented expected parameter type for `convergents`, if you're passing a tuple then it should be a pair of `int`s, but you're passing a pair consisting of an `int` and a `float`. I don't think that's supported. For SymPy, it would be helpful if you showed the relevant code. (How are you defining `sqrN`, and why are you using `//` rather than plain old `/`?) – Mark Dickinson Jul 13 '21 at 11:22
  • Also, are you aware that `math.sqrt(2140e225)` will still be a floating-point number, so will be rational (in fact, it'll even be an integer). That's not going to help you if you want to do exact arithmetic involving square roots; SymPy is likely the most reasonable way forward. – Mark Dickinson Jul 13 '21 at 11:46
  • Mark Dickinson, I updated question with info – Alexandr Dorofeev Jul 13 '21 at 14:58
  • If you want the entire continued fraction up to the part where it repeats, then for that size number you're asking for something that's infeasible to compute. You should expect the continued fraction for `1 / sqrt(N)`, for an arbitrarily chosen `N`, to be periodic with period of order of magnitude `sqrt(N)` (*very* roughly speaking). So that's going to be computable maybe up to `N = 10^16` or so. `2140e225` is _way_ beyond what's reasonable. – Mark Dickinson Jul 13 '21 at 16:46
  • To start with, you shouldn't calculate the CF of `1/√n` but `√n` and then just insert a 0 as the initial element by shifting all previously obtained coefficients to right. In CF world unshifting a zero means reciprocate. To obtain `√n` the best possible way, you may check [this answer of mine](https://stackoverflow.com/a/71306796/4543207) and also [this](https://medium.com/p/470a12044940) – Redu Aug 24 '22 at 15:21

1 Answers1

0

Just found one gist. It gives the same values as Wolfram's Alpha Convergents[]. However I'm not sure in precision after 5th element, because Wolfram only gives values up to that point.

import math
import itertools

def continued_fraction(n, d):
    while d:
        q, r = divmod(n, d)
        n, d = d, r
        yield q


def alternative_continued_fraction(n, d):
    gen = continued_fraction(n, d)
    p = next(gen)
    for q in gen:
        yield p
        p = q
    yield p - 1
    yield 1


def convergents(n, d):
    hh, kk, h, k = 0, 1, 1, 0
    for x in continued_fraction(n, d):
        hh, kk, h, k = h, k, h * x + hh, k * x + kk
        yield h, k

decimal.getcontext().prec = 500
print(list(convergents(1, Decimal(find_number).sqrt())))
Alexandr Dorofeev
  • 470
  • 1
  • 4
  • 12