prime numbers and fibonacci series

Question

Well let's say that we have two arrays in which we have a prime numbers sequence and a Fibonacci sequence e.g.

Fib = [1,1,2,3,5,8]

and

Pr = [2,3,5,7,11,13]

I need to find the first two numbers for which

(Fib[n1] % Pr [n1])==0 and (Fib[n2] % Pr [n2]) == 0

My code seems to be ok syntactically but is very slow when it comes to finding those numbers.Any help would be appreciated.

Code:

import math
def main():
    def Fib(n):
        v1, v2, v3 = 1, 1, 0   
        for rec in bin(n)[3:]:  
            calc = v2*v2
            v1, v2, v3 = v1*v1+calc, (v1+v3)*v2, calc+v3*v3
            if rec=='1':    v1, v2, v3 = v1+v2, v1, v2
        return v2
    def PrimeNum(n):
        P = ["",2]
        i = 3
        while len(P)!=(n+1):
            z = True
            if math.sqrt(i) == int(math.sqrt(i)):
                z = False
            else:
                j = 2
                while j<math.sqrt(i):
                    if i%j == 0:
                        z = False
                        break
                    j += 1
            if z == True:
                P.append(i)
            i += 2
        return P
    n1 = 0
    n2 = 0
    i = 1
    while Fib(i)%PrimeNum(i)[i] != 0:
        i+=1
    n1 , i = i , n1+1
    while Fib(i)%PrimeNum(i)[i] != 0:
        i+=1
    n2 = i
    print (n1+n2)
main()

There are tons of questions for this you could have searched for. — Pythonista, Dec 07 '16 at 03:33
I did search for this in other questions but i did not find any related post to this — ADR, Dec 07 '16 at 04:23
should be `while j<=math.sqrt(i):`, otherwise you include squares of primes, like 25. And your `PrimeNum(i)` recalculates all the `i-1` previous primes, instead of keeping the primes list and adding just one more prime to it (same as for the Fibonacci numbers). — Will Ness, Dec 07 '16 at 11:14
do not recalculate, *reuse*! do that and your code will *fly*. — Will Ness, Dec 07 '16 at 11:30
@WillNess j is a float number so by adding 1 i will never get the exact value of sqrt(5) so the equality doesnt affect the code i think but i'll check it out. — ADR, Dec 07 '16 at 12:48

eyllanesc · Answer 1 · 2016-12-07T04:29:41.017

import math


def Fib(n):
    v1, v2, v3 = 1, 1, 0
    for rec in bin(n)[3:]:
        calc = v2 * v2
        v1, v2, v3 = v1 * v1 + calc, (v1 + v3) * v2, calc + v3 * v3
        if rec == '1':    v1, v2, v3 = v1 + v2, v1, v2
    return v2


def isPrime(n):
    if n <= 1:
        return False
    elif n <= 3:
        return True
    elif (n % 2 == 0) or (n % 3 == 0):
        return False
    i = 5

    while i**2 <= n:
        if (n % i == 0) or ((n % (i+2)) == 0):
            return False
        i += 6
    return True


def PrimeNum(n):
    if n == 1:
        return 2
    else:
        i = 1
        counter = 1
        while counter < n:
            i += 2
            while not isPrime(i):
                i += 2
            counter += 1
        return i


def main():
    numbers = []
    n = 1
    while len(numbers) < 2:
        prime = PrimeNum(n)
        fib = Fib(n)
        if (fib % prime) == 0:
            numbers.append((n, prime, fib))
        n += 1
        print(n)

    print(numbers)

if __name__ == '__main__':
    main()

Output:

[(2160, 19009, 11582916825736584646975443653366322388921641167905180075281721851717008301903980437476277463379606054064279901980733340808520924624506359466585610824954684399084897922087504296705002387813173728581009860851320738456488209562894492182740289488332663607728717459344508234312101775586816401162288539598276213213951642360269055931032970100599429340189881639747365248788048662076913003130071131663108327512733069326012197090400460556913764055449808782413120), (3048, 27941, 4414103171422717935355962509499454349778622162400781459697552399266028998073052301861983592098091170325382922374854427709365699448595338643109175128231512749180137239630552350668736159301331619783913278988698914397585103122958242938481707043703775239751032229629764002510837997551552225537574210211945427154275570535341744487610518195389286251136162913276874919649834195343551175847369401914668124502580414188411921940167995410000017796803115148941554642484534655371013695818145603303564888547079314332213882853124811024475730089466430117462078329371553861038010496072275843988120946471902434970719656034482802077255302581089149251082976)]

yes it's a typo in the explanation of the algorithm not the algorithm itself i edited my first post sorry for that. — ADR, Dec 07 '16 at 04:26
@eyllanesc I tried running your code but it doesn't seem to work either it takes far too long to compute those numbers — ADR, Dec 07 '16 at 04:37
@ADR yes, the same problem of recalculating progressive prefixes of the primes sequence, instead of just adding them one by one, is still there. — Will Ness, Dec 08 '16 at 10:20

John Coleman · Accepted Answer · 2016-12-08T12:25:09.017

On Edit: As @WillNess pointed out, my original sieve-based method of generating primes was suboptimal. I have revised the prime generator.

With each and every n you are reconstructing all previous Fibonacci numbers in computing the nth Fibonacci number. This makes your algorithm for generating Fibonacci numbers at least quadratic in n -- not a good way to compute the terms of a linear recurrence.

Similarly with each and every n you are testing if n is prime by a crude trial division. It would make more sense to store already computed primes to stream-line that part, and only test odd numbers (after 2) for being the next prime.

Something like this:

from math import sqrt

def isPrime(k,known):
    """assumes k is odd and known includes odd primes <= sqrt(p)"""
    s = sqrt(k)
    for p in known:
        if k % p == 0:
            return False
        elif p > s:
            return True

def primes():
    yield 2
    yield 3
    yield 5
    known = [5] #known odd primes > 3. Iteration skips multiples of 3
    candidate = 7
    parity = 1

    while True:
        while not isPrime(candidate,known):
            candidate += (2 + 2*parity)
            parity = 1 - parity
        yield candidate
        known.append(candidate)
        candidate += (2 + 2*parity)
        parity = 1 - parity



def fib():
    yield 1
    a = 1
    b = 1
    while True:
        yield b
        a,b = b, a+b

def search(k, limit = 100000):
    """searches for the first k examples among first limit pairs"""
    hits = []
    for i,(f,p) in enumerate(zip(fib(),primes())):
        if i > limit:
            return "Not found"
        elif f % p == 0:
            hits.append((i,f,p))
            if len(hits) == k: return hits

When search(2) is invoked, it almost instantly returns:

[(2159, 11582916825736584646975443653366322388921641167905180075281721851717008301903980437476277463379606054064279901980733340808520924624506359466585610824954684399084897922087504296705002387813173728581009860851320738456488209562894492182740289488332663607728717459344508234312101775586816401162288539598276213213951642360269055931032970100599429340189881639747365248788048662076913003130071131663108327512733069326012197090400460556913764055449808782413120, 19009), (3047, 4414103171422717935355962509499454349778622162400781459697552399266028998073052301861983592098091170325382922374854427709365699448595338643109175128231512749180137239630552350668736159301331619783913278988698914397585103122958242938481707043703775239751032229629764002510837997551552225537574210211945427154275570535341744487610518195389286251136162913276874919649834195343551175847369401914668124502580414188411921940167995410000017796803115148941554642484534655371013695818145603303564888547079314332213882853124811024475730089466430117462078329371553861038010496072275843988120946471902434970719656034482802077255302581089149251082976, 27941)]

For fun:

>>> search(3)
[(2159, 11582916825736584646975443653366322388921641167905180075281721851717008301903980437476277463379606054064279901980733340808520924624506359466585610824954684399084897922087504296705002387813173728581009860851320738456488209562894492182740289488332663607728717459344508234312101775586816401162288539598276213213951642360269055931032970100599429340189881639747365248788048662076913003130071131663108327512733069326012197090400460556913764055449808782413120, 19009), (3047, 4414103171422717935355962509499454349778622162400781459697552399266028998073052301861983592098091170325382922374854427709365699448595338643109175128231512749180137239630552350668736159301331619783913278988698914397585103122958242938481707043703775239751032229629764002510837997551552225537574210211945427154275570535341744487610518195389286251136162913276874919649834195343551175847369401914668124502580414188411921940167995410000017796803115148941554642484534655371013695818145603303564888547079314332213882853124811024475730089466430117462078329371553861038010496072275843988120946471902434970719656034482802077255302581089149251082976, 27941), (27093, 9154611756214756724173885190408003872449591269073102605228456275421208215932565362931929654546649304023236344043061900343945422952312397120822618708723071089543994689417534274219382267837418988544430955958634578186146315678653949605935409450055790001104946147407316311473901552255789990255560458743971937801679393705446136261518016083185450683814066644274897534221110708848082023454498620988462247604771107063854440837945572552841406526401192557118031116749129168215879638214765494261573532277267248651864394801267996413906424450730419095713253333326504619013654894156288241737988067918591507737206648058569934516619420144623505377301435668514418593118878788775849264442773570236710982356483749339325612961955144121852617030353522266206195503391667491427262341962604316852915673071607898002570374470506997388337604296115573269049069634099324411619052448132528877435091667009229816989818556680426118774810238582636027323802599258271991479814852819141878003265064480606587244021927725713293228979923827422432692586158272135063549520191382577769000483431081463036690550524564226586425639020224650313666850943741943605246548446573870629218059701374187171356525816670387118159122398327173724524855658200110122522835176190387328244551014672164462460351372690457712069935724023091919036138410906797524707392410554446350077959242405947099126465838760766371777403033146251749502951971723164017865200568964463774297505017369007733498550947322388608760212335940368621783904555284384376776656261924863794343796815840394007873484565022848579462774052069439946777630573102979496590395316870340899551813282502986710070529378861792991385680534480763081222794445860234260101963233469819968478267781058165962215892733089244018231780105097548621006518924522684510451706606468473417782105562726421667458599325050626566145642109892327168775105402123881849220660482291948715812488950704140069539809853344947732388220095254658324470691606583872184849553953678352230667797760550750875445274878011552908141893541049214761478631614429073173072195839699424914223620432036633802986310860934484009114489007981128832162565709338385636160750356924485343604148697926502326166202557944994913633051828299661050405190916615302885951550480867328010617078304297055489653587862296405779015531981276461702240597546551137530304640351560899721596431710387198320317078019194398346106906325232766448980878644774607110433151004714381896844981823117568831262694359188988036994358418407531621088012313911112990787883908629473366729105235570919108365315772383671485216041535895199802233506396125361100590961904686232499367884224371237416713503169202893109931633850881704747136301856000860027043777582090218562843855449359939394306400414888273691232591121594524381819354582423430557415390250709297250916474909569851542211629963074881157328576497818802346387745923866756229123554832128989541520169936300953836300235389379093109442248393586494022605971670012623627311681211479804894428255735216832312332374587671329106796684022978470999078866754294893462641324866080667457196130998327268153752398794195850817076787793178090421859530598534071035337698763058968997892119530295325366370181024862996688996982218234484156138212120039174745114551380799834389352242746434596156561884802858552848039130654598584670423370789192360682663573402610067394573907023581267337444561887066225448640531109903737583527982043515262559110647745958317794321134493638412203902264050559048131283893050978643459005295587397025818844745289630423132369876077023573075795101695741941448957492692770541027257813143549818083247077466563300619633483770879467622427653718311051255568315886707370049330370551248083587163093648392436402538080835372938539661876694968429300607645700887451607717345890164088866738067437705736191409017874160533226574366616228364590208669132008603977373406173169617176940813788990290600726024523687722649543950416403252632629816054630374545316484236530285427928212536677046788146607820784957497688177505215507738764368094945972293925768102942894343922075951950584747843163041482179651582945321367224573871167076773239474546718093054690180737239167052175372000610893947129090802650539145357698139189391292101318833102423529064084510914322440126808385923038318201245755676704650611731596292583717304379370256996717858556466693113975077196901855204797418153369689670038673829781023407381326668723988156703032126856483599703246194765161592028779525610900396244517331650235058027054924312880572459575090013084993690584284681466416902902178987116510671639573471785558879473018150310384889626008363718980810299980477673370182512022494547611737471877682888208577935567297869927069037836778041611528430972865356012323596536257940865191759958191141927469705717229297160293645616965252935768838221069724045314321468122778102190272203721595427118234293979687595632013680647352722589129686173006256611741794504701897170430155303809328057620955332830967596645368594105246945340269177608454700950981542453098799391204902996079682585726077788735323172833669157208933105328320688115165403826545598927201627249410032419171147416036060482494310313973399963662345404911693098014230729136599941932163456611595895059104775163455217728743575034276058374708178334997147073744767691582680892897199796122494961096167960722019157173361146936811991971390168385637304432676953904476762445279910913171572513177072541852580653529021831314795629186219681562085784871790357946471666771284027378207178520108478298275531913586518715536383856596319509987179659668437785231444827668882063238119727411375182795461735962704177295871116919234266198086260511993082345239690462198949212239357280276703513750501938682329025159521973018772292717925652964949167, 313721)]

When I increase limit, I am able to find a 4th example, where the prime in question is 3523969. The sequence of indices of such pairs is A075702 at the Online Encyclopedia of Integer Sequences. This reference doesn't indicate whether or not the sequence is known to be infinite.

Thank you very much your code works flawlessly in a fraction of the time that mine did. I'll just try to understand it right now through your thorough explanation :P — ADR, Dec 07 '16 at 05:06
@ADR The second version I added is a little less cluttered and runs slightly faster. You can always mark an answer as "accepted" if you want. — John Coleman, Dec 07 '16 at 05:15
when I try to run it with : searchTo(5000) it gives me "not found" as a result but in the previous version I got the results I wanted, why is that happening? — ADR, Dec 07 '16 at 05:27
The primes in question are > 5000. The parameter controls what range of numbers I am sieving when I generate the primes. Note that in the sample code I used `searchTo(50000)` rather than `searchTo(5000)`. Of course, before you know the answer there would be a certain amount of trial and error, though even something like `searchTo(100000)` isn't too bad. The sieving part of the code could still use some optimization. — John Coleman, Dec 07 '16 at 05:31
but the second prime for which F[n]%Pr[n] == 0 is the 3047th in the prime's number series so searchTo(5000) should work but it actually doesn't. I think I got it . It doesnt search for the position of the prime number but for the number itself as it is the number : 27941 — ADR, Dec 07 '16 at 05:38
The 3047th prime is 27941 -- n is a bound on the *size* of the prime -- not its index. Its purpose is to tell the generator how much internal state it has to maintain. `n` has a different meaning in the revised code than in the original code, which does make it somewhat confusing. — John Coleman, Dec 07 '16 at 05:57
oh, no, the "generator sieve" is a suboptimal trial division inferior to the OP's "crude trial division", which is in fact not crude, but optimal! The problem there was the recalculation of primes, just as you've noted for the Fibonacci numbers. --- BTW it'd only be quadratic if fib(n) was O(1) but it is O(log(n)); so it's even worse. — Will Ness, Dec 07 '16 at 11:18
as for primes, your code is O(n^2), in `n` primes produced, while the optimal trial division is O(n^1.5) and true sieve is O(n log(n)); the OP code, because of the recalculations, was about O(n^2.5). So primes actually might dominate the fibonaccis O(n^2 log(n)) there -- except for the bignum overhead. — Will Ness, Dec 07 '16 at 11:33
@WillNess Thanks for the heads-up. In my original approach I just used an implementation of a crude sieve that I had previously written when playing around with the Sieve of Eratosthenes. Based on your comment, I scrapped that generator in favor of a new one. It is now able to find the first 2 examples almost instantly, and can find the first 3 in about the same time my first approach took to find 2. Thanks again. — John Coleman, Dec 07 '16 at 12:19
glad to be of help. I misspoke, the OP trial division wasn't optimal, it used all numbers and not just primes, but still it stopped at the sqrt which is more important (it's basically sqrt vs log (for using only primes instead of all numbers, or odds)). — Will Ness, Dec 07 '16 at 14:10
there's still one major redundancy in your updated `primes()`. Would it interest you to find out what it is? :) — Will Ness, Dec 08 '16 at 10:18
@WillNess Well, my way of alternating between adding 2 and adding 4 is a bit silly since I "save" an evaluation of `candidate % 3` at the cost of an evaluation of `candidate % 6` -- there are much more efficient ways to alternate. But -- that is a coding infelicity rather than a major redundancy. I am of course interested in improving it and I'll think of it a bit more, but am about to plunge into finals week where I work so I won't have much time in the next 10 days or so. Number theory isn't my area, so I am aware that code is on the naïve side. — John Coleman, Dec 08 '16 at 12:18
I'll give you a hint: it's about its space use, not time. I won't spoil your fun. (and yes, it implements optimal trial division - not the sieve of Eratosthenes - but that's besides the point). — Will Ness, Dec 08 '16 at 12:27
@WillNess I could boot-strap. Have a "fast" generator which is the main generator. It could maintain an internal "slow" generator which produces the primes to check against. The internal generator need only be invoked when the fact generator passes another perfect square. Total primes stored will only be about the square root of what I am currently doing. It does involve a space vs. time trade-off but wont increase big-O. — John Coleman, Dec 08 '16 at 13:28
you got it. :) and it could actually run faster too. in [another, related setting](http://stackoverflow.com/questions/2211990/how-to-implement-an-efficient-infinite-generator-of-prime-numbers-in-python/10733621#10733621), it did. — Will Ness, Dec 08 '16 at 17:43

prime numbers and fibonacci series

2 Answers2