一般数体ふるい法で使うふるいの実装トライアル

一般数体ふるい法に使うふるいのプログラムを作ってみました。ただ、難しくてどう実装すればいいのかわからない計算式もあったので、できる範囲だけでの実装です。4桁くらいの小さい数までしかしっかりと動かせない感じです。

import numpy as np
import math
import random



def ugojod(a,b):
    if a<0:
        a*=-1
    if b<0:
        b*=-1
    if a>b:
        a,b=b,a
    while a>0:
        c=int(b)//int(a)
        a,b=b-c*a,a
        #print(a,b)
    return b
def kaihei(a,sisuu):
    #print(a,'a')
    ns=[]
    while a>0:
        ns+=[a%10**sisuu]
        a//=10**sisuu
    #print(ns)
    m=len(ns)
    tm=0
    c=0
    ans=0
    for i in range(m):
        b=0
        #print(c,ans,tm)
        c=c*10**sisuu+ns[m-i-1]
        for j in range(10):
            #print(j,tm,b)
            tm+=1
            b+=1
            if b*tm<=c:
                1
            else:
                tm-=1
                b-=1
                break
            
        #print(b)
        ans=ans*10+b
        c-=b*tm
        tm=(tm+b)*10
    #print(c,ans,'ans')
    return c,ans 
def xgcd(a,b):
    x0,y0,x1,y1=1,0,0,1
    while b!=0:
        q,a,b=a//b,b,a%b
        x0,x1=x1,x0-q*x1
        y0,y1=y1,y0-q*y1
    return a,x0,y0
def modinv2(a,m):
    g,x,y=xgcd(a,m)
    if g!=1:
        return 'not exist'
    else:
        return x%m 
def bezulis(baselist,amarilist):
    b=baselist
    a=amarilist
    sh=a[0]
    sab=b[0]
    for i in range(1,len(b)):
        go=ugojod(sab,b[i])
        if go>1:
            if ugojod(a[i]-sh%b[i],go)==go:
                sh+=sab*(a[i]-sh%b[i])*modinv2(sab//go,b[i]//go)//go
                sab*=b[i]//go
            else:
                continue
        else:
            sh+=sab*(a[i]-sh%b[i])*modinv2(sab,b[i])
            sab*=b[i]
        #print(sh,sab)
    return sh,sab
    
def bekiyolists(bunkaiN,tei0,bekilist):
    #bunkaiN=497
    #tei=4
    #beki=12
    beki=len(bekilist)
    bekiyolist=[0]*beki
    blist=[]
    tei=tei0
    l=0
    bekiloglist=[[]]*beki
    mx=0
    #print(bekilist)
    for j in bekilist:
        #print(j)
        beki2=j
        #if l%100==1:
            #print(l)
        i=0
        listes=[]
        while beki2>0:
            listes.append(beki2%2)
            #blist.append(beki2%2)
            beki2=beki2//2
            i+=1
        bekiloglist[l]=listes
        if mx<i:
            mx=i
        l+=1
    #print(bekiloglist)
    result=[1 ]*beki
    for i in range(mx):
        for j in range(beki):
            if len(bekiloglist[j])>i:
                #print(result[j],tei,bekiloglist[j][i],bunkaiN)
                result[j]=(result[j]*tei**bekiloglist[j][i])%bunkaiN
        #tei=bigprod(tei,tei)%bunkaiN
        tei=(tei**2)%bunkaiN
        #print(tei,'tei',bunkaiN)
    bekiyolist=result
    #print(result)
    #print(beki,blist)
#blist.reverse()
        #result=1

        #for bv in blist :
            #if bv==1:
                #result=(result*tei)%bunkaiN
            #tei=(tei*tei)%bunkaiN
        #bekiyolist[k]=result
        #k+=1

        #print(tei,result)
    return bekiyolist 
def fermat(a):
    c=3
    if a==2 or a==3:
        return True
    if bekiyolists(a,2,[(a-1)])[0]==1 and bekiyolists(a,3,[(a-1)])[0]==1:
        return True
    return False
def millerrabin(a):
    if a==2:
        return True
    if a>2 and a&1==0:
        return False
    s,t=0,a-1
    while t&1==0:
        s,t=s+1,t>>1
    b=[2,7,61]
    ret=False
    for i in b:
        if i>a:
            break
        if bekiyolists(a,i,[t])[0]==1:
            #print(i,s,t,a,ret)
            return True
            ret= True
        tt=t
        for j in range(0,s):
            if bekiyolists(a,i,[tt])[0]==a-1:
                #print(i,s,j,tt,a,ret)
                return True
                ret=True
            tt*=2
            
    return ret 
def heihoufind2(a,bl,pnl=[]):
    pro=1
    for i in range(len(bl)-3):
        pro*=bl[i]
        if pro>a:
            bl=bl[0:i+3]
            break
    cl=[]
    dl=[bekiyolists(bl[i],a,[(bl[i]-1)//2])[0] for i in range(len(bl))]
    #print([bekiyolists(bl[i],a,[(bl[i]-1)//2])[0] for i in range(len(bl))],bl)
    if len(pnl)>0:
        for i in range(len(bl)):
            if dl[i]==1:
                cl+=[pnl[i][0]]
            elif dl[i]+1==bl[i]:
                cl+=[pnl[i][1]]
            else:
                cl+=[[]]
    else:
        for i in range(len(bl)):
            ccl=[]
            for j in range(1,bl[i]):
                if bekiyolists(bl[i],j,[(bl[i]-1)//2])[0]==dl[i]:
                    ccl+=[j]
            cl+=[ccl]
    #print(cl,bl,dl,'dl')
    mo=cl[0]
    l=len(bl)
    pro=bl[0]*bl[1%l]*bl[2%l]
    d=[i for i in range(pro)if i%bl[0] in cl[0] and i%bl[1] in cl[1] and i%bl[2] in cl[2]]#+[3319]
    e=[2 for i in range(len(d))]#+[2]
    #print(d,e)
    #print(8 in [2,4,8,10])
    ret=[]
    #print([bekiyolists(bl[j],274,[(bl[j]-1)//2])[0] for j in range(l)],274)
    for k in range(8):
        pro=bl[k%l]*bl[(k+1)%l]*bl[(k+2)%l]
        dd=[]
        ee=[]
        fl=0
        for i in range(len(d)):
            if (k+2-e[i])%l==0 and d[i]%bl[(k+2)%l] in cl[(k+2)%l]:
                #print(d[i],[bekiyolists(bl[j],d[i],[(bl[j]-1)//2])[0] for j in range(l)],dl)
                if dl==[bekiyolists(bl[j],d[i],[(bl[j]-1)//2])[0] for j in range(l)]:
                    ret+=[d[i]]
                    fl=0
                    #print(xs)
            if d[i]%bl[(k+2)%l] in cl[(k+2)%l]:
                dd+=[d[i]]
                ee+=[(e[i])%l]
                pp=d[i]//pro
            f=[pp*pro+(d[i]+j*bl[(k)%l]*bl[(k+1)%l])%pro for j in range(bl[(k+2)%l])]
            f=list(set(f))
            f=[f[j] for j in range(len(f))if f[j]%bl[(k+2)%l] in cl[(k+2)%l] and f[j]!=d[i] and f[j]<=a]
            dd+=f
            ee+=[(k+3)%l]*len(f)
        d=dd
        e=ee
        f=list(set(d))
        gg=[]
        for i in range(len(f)):
            g=[e[j] for j in range(len(e))if d[j]==f[i]]
            ggg=[(k+2-g[j])%l for j in range(len(g))]
            gg+=[g[ggg.index(max(ggg))]]
        d=f
        e=gg
        #print(dd,ee,k,bl[(k+2)%l],(k-2)%l)
        #print([d[j] for j in range(len(e))if e[j]==2],bl[(k+2)%l])
        #print('lend',len(dd))
        if fl==1:
            break
        if k%4==3:
            dd=[]
            ee=[]
            fl=0
            for i in range(len(d)):
                if d[i] in cl[(k+3)%l] or (k+2-e[i])%l>5:
                    dd+=[d[i]]
                    ee+=[e[i]]
            d=dd
            e=ee
        #print('lend',len(dd))
    #print(ret,l)
    #print(mo)
    ret2=[]
    if len(ret)==0:
        return a,1
    for i in range(len(ret)):
        if a%ret[i]==0 and a>ret[i]:
            if kaihei(a//ret[i],2)[0]!=0:
                continue
            ret2+=[ret[i]]
    print(ret2)
    if len(ret2)==0:
        return a,1
    ans=max(ret2)
    return ans,a//ans 
def funcchoose(a):
    b=kaihei(a,2)[1]
    c=[a//b**2]
    c+=[(a-c[0]*b**2)//b]
    c+=[a-c[0]*b**2-c[1]*b]
    print(b,c)
    if c[1]%2==1:
        c[1]-=1
        c[2]+=b
    cc=[]
    dx=0
    if c[1]%2==0:
        cc+=[1,0,c[2]-(c[1]//2)**2]
        dx=b-c[1]//2
    else:
        cc+=[4,0,4*c[2]-(c[1])**2]
        dx=2*b-c[1]
    print(cc,dx)
    return cc,dx
def quadquad(a,ls,lim=1000):
    b=kaihei(a,2)[1]
    limit=lim
    if b<lim:
        limit=b
    e=[]
    f=[]
    for i in range(limit):
        c=a-(b-i)**2
        d=heihoufind2(c,ls)
        e+=[d[0]]
        f+=[d[1]]
        #print(c,d)
        #print(xs)
        if i%100==99:
            print(i)
    g=e.index(min(e))
    return [1,0,e[g]],b-g,kaihei(f[g],2)[1] 
def soide2(xl,n):
    idels=[]
    alpls=[]
    prils=[]
    irange1=1
    irange2=200
    jrange=200
    krange=500
    prirange=80000
    pri2range=2000
    nidels=[]
    nprils=[]
    count=0
    for i in range(irange1,irange2):
        if i%10==9:
            print('i',i)
        for j in range(1,jrange):
            a=xl[0]*i**2+xl[2]*j**2
            if i==1 and j==1:
                print(a,i,j,xl,fermat(a),millerrabin(a),([a] in prils)==False,prils)
            if 1:
                if 1:
                    bb=1
                    if fermat(a) and millerrabin(a) and a<prirange and( a in prils)==False:
                        alpls+=[bb]
                        #idels+=[[i,j,1,0]]
                        idels+=[[i,j,0,1,0]]
                        prils+=[a]
                        #count+=1
                        #if bekiyolists(a,bb,[a-1])[0]!=1:
#                            continue
                    kkrange=krange
                    if kkrange>a:
                        kkrange=a
                    l=0
                    for l in range(1,kkrange):
                        if (-xl[2]*modinv2(xl[0],a))%a==l**2%a:#check
                            ll=l
                            #count+=1
                            break
                    else:
                        continue
                    for k in range(1,krange):
                            kk=(-k*ll)+a#kk**2+k/ll**2==0
                            if (kk**2*xl[0]+k**2*xl[2])%a==0:
                                #count+=1
                                p=(kk**2*xl[0]+k**2*xl[2])//a
                                #print(p)
                        #print(i,j,k,a)
                                if p>=2 and fermat(p) and millerrabin(p) and p<prirange and( p in prils)==False:
                                    bbb=1
                                    alpls+=[bbb]
                                    idels+=[[kk,k,0,i,j]]
                                    prils+=[p]
                                    count+=1
    print(count,'count')
    for i in range(2,pri2range):
        if fermat(i) and millerrabin(i) and i not in prils:
            alpls+=[0]
            idels+=[[(i*modinv2(1,n))%n+xl[0],0,xl[0],1,0]]
            prils+=[i]
    print(count,'count')
    print(idels[0:20],alpls[0:20],prils[0:20])
    print(len(idels),len(alpls),len(prils))
    print(len(idels))
    print(list(sorted([abs(prils[k]) for k in range(len(prils))]))[0:90])
    return idels,alpls,prils 
def idesieve(ide,fx,n,xx,yy=1):
    ils=ide[0]#[0:20]
    als=ide[1]#[0:20]
    pls=ide[2]#[0:20]
    irange1=1
    irange2=150
    jrange1=2
    jrange2=500
    prirange=3000
    #print(pls)
    #print(xs)
    spls=[]
    anspri=[]
    anspls=[]
    anspair=[]
    count=0
    pri=[s for s in range(2,prirange)if fermat(s) and millerrabin(s)]
    for i in range(irange1,irange2):
        
        if i==0:
            continue
        for j in range(jrange1,jrange2):
            if j==0:
                continue
            if ugojod(i,j)>1:
                continue
            c=[]
            #a=i**2+j**2
            a=fx[0]*i**2+fx[2]*j**2
            g=[]
            alp=[]
            for k in range(len(pls)):
                if ils[k][2]==0 and ils[k][4]==0 and a%abs(pls[k])==0 :#and (i+j*als[k])%abs(pls[k])==0:
                #if a%pls[k]==0:
                    c+=[ils[k]]
                    #d*=fx[0]*ils[k][0]*yy+fx[2]*ils[k][1]*xx
                    g+=[pls[k]]
                    alp+=[als[k]]
                    continue
                if ils[k][2]!=0 and ils[k][4]==0 and a%abs(pls[k])==0:#and (i+j*als[k])%abs(pls[k])==0:
                #if a%pls[k]==0:
                    c+=[ils[k]]
                    #d*=fx[0]*ils[k][0]*yy+fx[2]*ils[k][1]*xx
                    g+=[pls[k]]
                    alp+=[als[k]]
                    continue
                
                if ils[k][4]!=0 and a%abs(pls[k])==0:
                #if a%pls[k]==0:
                    
                        
                    c+=[ils[k]]
                    g+=[pls[k]]
                    alp+=[als[k]]
            if len(c)>0:
                s=i*yy+j*xx
                ss=s
                prilist=[]
                plslist=[]
                ilslist=[]
                pplslist=[]
                alplist=[]
                tatamils=[1]
                d=1
                rd=1
                for k in range(len(pri)):
                    while s%pri[k]==0 and s!=0:
                        s//=pri[k]
                        prilist+=[pri[k]]
                aa=a
                for k in range(len(g)):
                    if g[k]>0 and c[k][2]==0 and c[k][4]==0:
                        while (aa)%abs(g[k])==0 and aa!=0 and abs(g[k])!=yy:
                                d*=c[k][0]*yy+c[k][1]*xx
                                ss*=yy
                                prilist+=[yy]
                                plslist+=[c[k][0]*yy+c[k][1]*xx]#[g[k]]
                                #aa*=yy
                                aa//=g[k]
                        #print(g[k])
                            #plslist+=[nor]#[g[k]]
                            #prilist+=[nor2]
                                ilslist+=[c[k]]
                                pplslist+=[g[k]]
                                #alplist+=[alp[k]]
                    if g[k]>0 and c[k][2]==0 and c[k][4]!=0:
                        while (aa)%abs(g[k])==0 and aa!=0 and abs(g[k])!=yy:
                                #d*=c[k][0]*c[k][3]*yy**2-c[k][1]*c[k][4]*xx**2+(-c[k][0]*c[k][4]+c[k][1]*c[k][3]*xx*yy)
                                d*=c[k][0]*yy+c[k][1]*xx
                                ss*=c[k][3]*yy+c[k][4]*xx
                                #ss*=yy
                                #ss*=c[k][3]**2*yy**2-c[k][4]**2*xx**2
                                prilist+=[c[k][3]*yy+c[k][4]*xx]
                                plslist+=[c[k][0]*yy+c[k][1]*xx]#[g[k]]
                                aa//=g[k]
                        #print(g[k])
                            #plslist+=[nor]#[g[k]]
                            #prilist+=[nor2]
                                ilslist+=[c[k]]
                                pplslist+=[g[k]]
                                #alplist+=[alp[k]]
                    
                if i%10==9and j==3 :#and (d-ss*rd)%n==0:
                    print(i,j,s,aa,d,ss,(d*aa-ss*s)%n,prilist,plslist,pplslist)
                #if s==1 and aa==1 and (d-ss)%n==0 and d!=ss:
                if (d*aa-ss*s)%n==0 and d*aa!=ss*s:
                    #print(i+j*x,i,j,x,c,d,prilist,aa)
                    print(d,ss,i,j,ss,a,s,aa,'ch',plslist,prilist,pplslist,ilslist)
                    anspri+=[prilist]
                    anspls+=[plslist]
                    anspair+=[[d*aa,ss*s]]
            #print(c)
    print(len(anspair))
    print(len(anspri),len(anspls),len(anspair))
    print('count',count)
    return anspair,anspri,anspls
                
a=131*1031#1565912117761
b=3
c=3
d=35
bl=[i for i in range(2,500)if fermat(i) and millerrabin(i)]

qx=quadquad(a,bl,50)
print(a)
#print(fx)
print(qx)
print(kaihei(a,2))
so=soide2(qx[0],a)#qx[0]
#print(so)
#print(xs)
idesieve(so,qx[0],a,qx[1],qx[2])

苦戦の跡がたくさんありますが、何故苦戦しているかと言えば、素イデアルのノルムが小さい素イデアルを作るのが、分解したい数が大きくなるにつれて難しくなるからです。。