3890 | | for i in range(ceil(((2*H_q+1)**r)/2)): |
3891 | | koeffs = Z(i).digits(base=2*H_q+1) |
3892 | | koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s |
3893 | | P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0)) |
3894 | | for Q in tors_points: # P + torsion points (includes 0) |
3895 | | tmp = P + Q |
3896 | | if not tmp.is_zero(): |
3897 | | x = tmp[0] |
3898 | | if x.is_integral(): |
3899 | | xs.add(x) |
| 3892 | N=H_q |
| 3893 | |
| 3894 | def use(P): |
| 3895 | """ |
| 3896 | Helper function to record x-coord of a point if integral. |
| 3897 | """ |
| 3898 | if not P.is_zero(): |
| 3899 | xP = P[0] |
| 3900 | if xP.is_integral(): |
| 3901 | xs.add(xP) |
| 3902 | |
| 3903 | def use_t(R): |
| 3904 | """ |
| 3905 | Helper function to record x-coords of a point +torsion if integral. |
| 3906 | """ |
| 3907 | for T in tors_points: |
| 3908 | use(R+T) |
| 3909 | |
| 3910 | # We use a naive method when the number of possibilities is small: |
| 3911 | |
| 3912 | if r==1 or N<=10: |
| 3913 | for P in multiples(mw_base[0],N+1): |
| 3914 | use_t(P) |
| 3915 | return xs |
| 3916 | |
| 3917 | # Otherwise is is very very much faster to first compute |
| 3918 | # the linear combinations over RR, and only compte them as |
| 3919 | # rational points if they are approximately integral |
| 3920 | |
| 3921 | def is_approx_integral(P): |
| 3922 | return (abs(P[0]-P[0].round()))<0.001 and (abs(P[1]-P[1].round()))<0.001 |
| 3923 | |
| 3924 | RR = RealField() #(100) |
| 3925 | ER = self.change_ring(RR) |
| 3926 | Rgens = [ER.lift_x(P[0]) for P in mw_base] |
| 3927 | for i in range(r): |
| 3928 | if abs(Rgens[i][1]-mw_base[i][1])>abs((-Rgens[i])[1]-mw_base[i][1]): |
| 3929 | Rgens[i] = -Rgens[i] |
| 3930 | for ni in cartesian_product_iterator([range(-N,N+1) for i in range(r-1)]+[range(N+1)]): |
| 3931 | RP=sum([ni[i]*Rgens[i] for i in range(r)],ER(0)) |
| 3932 | for T in tors_points: |
| 3933 | if is_approx_integral(RP+ER(T)): |
| 3934 | P = sum([ni[i]*mw_base[i] for i in range(r)],T) |
| 3935 | use(P) |
| 3937 | |
| 3938 | # Below this point we keep earlier experimental code which |
| 3939 | # is slower when either N or r is large |
| 3940 | |
| 3941 | # if r==2: |
| 3942 | # for P in multiples(mw_base[0],N+1): |
| 3943 | # for Q in multiples(mw_base[1],N+1): |
| 3944 | # for R in set([Q+P,Q-P]): |
| 3945 | # use_t(R) |
| 3946 | # return xs |
| 3947 | # |
| 3948 | # if r==3: |
| 3949 | # for P in multiples(mw_base[0],N+1): |
| 3950 | # for Q in multiples(mw_base[1],N+1): |
| 3951 | # PQ = [P+Q,P-Q] |
| 3952 | # PQ = set(PQ + [-R for R in PQ]) # {\pm P\pm Q} |
| 3953 | # for R in multiples(mw_base[2],N+1): |
| 3954 | # for S in PQ: |
| 3955 | # use_t(R+S) |
| 3956 | # return xs |
| 3957 | # |
| 3958 | # # general rank |
| 3959 | # mults=[list(multiples(mw_base[i],N+1)) for i in range(r)] |
| 3960 | # for i in range(r-1): |
| 3961 | # mults[i] = [-P for P in mults[i] if not P.is_zero()] + mults[i] |
| 3962 | # for Pi in cartesian_product_iterator(mults): |
| 3963 | # use_t(sum(Pi,self(0))) |
| 3964 | # return xs |
| 3965 | # |
| 3966 | # older, even slower code: |
| 3967 | # |
| 3968 | # for i in range(ceil(((2*H_q+1)**r)/2)): |
| 3969 | # koeffs = Z(i).digits(base=2*H_q+1) |
| 3970 | # koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s |
| 3971 | # P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0)) |
| 3972 | # for Q in tors_points: # P + torsion points (includes 0) |
| 3973 | # tmp = P + Q |
| 3974 | # if not tmp.is_zero(): |
| 3975 | # x = tmp[0] |
| 3976 | # if x.is_integral(): |
| 3977 | # xs.add(x) |
| 3978 | # return xs |