check_disc for mod-jon

%cython from sage.rings.integer cimport Integer from libc.stdint cimport int32_t, uint32_t, int64_t, uint64_t ctypedef int32_t int_cur[10] cdef mpz_t s,t,v def test(E): cdef int_cur realE N = ZZ(79)#ZZ(E.conductor().norm()) E = E.ainvs() cdef int32_t i for i in range(5): realE[2*i] = E[i][0] realE[2*i+1] = E[i][1] pri = N.prime_divisors() cdef uint32_t factLen = len(pri) cdef uint32_t *fact,*factEnd fact = <uint32_t *>sage_malloc(sizeof(uint32_t)*factLen) for i in range(factLen): fact[i] = pri[i] factEnd = fact+factLen mpz_init(t) mpz_init(s) mpz_init(v) cdef int32_t k = 10000000 w = walltime() for i in range(k): ret = check_disc(realE,N, fact, factEnd) print walltime(w)/k*10**9,'ns' mpz_clear(t) mpz_clear(s) mpz_clear(v) sage_free(fact) return ret cdef bint check_disc(int_cur ainv, uint32_t N, uint32_t *fact, uint32_t *factEnd): cdef int32_t b21, b22, b41, b42, b61, b62, b81, b82 cdef int64_t temp b21 = ainv[0]+ainv[1] b22 = ainv[1]*(ainv[0]+b21)+(ainv[3]<<2) b21 += ainv[2]<<2 b81 = b22*ainv[9] b82 = b21*ainv[9]+b22*ainv[8]+b81 b81 += b21*ainv[8] b41 = ainv[1]*ainv[5] b42 = b41+ainv[0]*ainv[5]+ainv[1]*ainv[4] b41 += ainv[0]*ainv[4] temp = ainv[7]*b42 b82 -= temp + ainv[6]*b42 + ainv[7]*b41 b81 -= temp + ainv[6]*b41 b41 += ainv[6]<<1 b42 += ainv[7]<<1 b61 = ainv[4]+ainv[5] b62 = ainv[5]*(ainv[4]+b61) temp = ainv[3]*b62 b82 += temp + ainv[2]*b62 + ainv[3]*b61 b81 += temp + ainv[2]*b61 b61 += ainv[8]<<2 b62 += ainv[9]<<2 temp = ainv[7]*ainv[7] b81 -= temp + ainv[6]*ainv[6] b82 -= temp + ((ainv[6]*ainv[7])<<1) cdef int64_t tmp1,tmp2 tmp1 = b22*b42 tmp2 = tmp1+b22*b41+b21*b42-3*b62 tmp1 += b21*b41-3*b61 temp = b62*tmp2 tmp2 = 9*(temp+b62*tmp1+b61*tmp2) tmp1 = 9*(temp+tmp1*b61) mpz_set_si(s,tmp1) mpz_set_si(t,tmp2) tmp1 = b41*b41 tmp2 = b42*b42 temp = (b42+3*b41)*tmp2 tmp2 = temp+(tmp2+3*tmp1)*b42<<3 tmp1 = temp+tmp1*b41<<3 if tmp1 >= 0: mpz_sub_ui(s,s,tmp1) else: mpz_add_ui(s,s,-tmp1) if tmp2 >= 0: mpz_sub_ui(t,t,tmp2) else: mpz_add_ui(t,t,-tmp2) temp = b21*b21+b22*b22 b22 = ((b21<<1)+b22)*b22 b21 = temp tmp1 = b81 tmp2 = b82 temp = b22*tmp2 tmp2 = temp+tmp1*b22+tmp2*b21 tmp1 = temp+b21*tmp1 if tmp1 >= 0: mpz_sub_ui(s,s,tmp1) else: mpz_add_ui(s,s,-tmp1) if tmp2 >= 0: mpz_sub_ui(t,t,tmp2) else: mpz_add_ui(t,t,-tmp2) mpz_sub(v,s,t) mpz_mul(v,v,t) mpz_addmul(v,s,s) if mpz_mod_ui(t,v,N): return 0 while fact != factEnd: while not mpz_mod_ui(t,v,fact[0]): mpz_divexact_ui(v,v,fact[0]) fact += 1 if mpz_sgn(v) < 0: mpz_neg(v,v) return mpz_root(v, v, 12ul) 

K.<a> = NumberField(x**2-x-1) #R.<x> = K[] 

E = EllipticCurve(K,[1,0,0,-784,-8515]) 

bb=E.b_invariants() 

E.discriminant()+bb[0]**2*bb[3] 

-623024

-623024

(-1568)^3*8 

-30840979456

-30840979456

147/21 

7

7

E.discriminant() 

147

147

E.conductor() 

Fractional ideal (21)

Fractional ideal (21)

(-1568+3*34060)*9*(-34060) 

-30841602480

-30841602480

test(E) 

186.378884315 ns
False

186.378884315 ns
False

 

x = 5*a+2 y = 7*a-3 

len(ZZ(E.discriminant()[1]).bits()) 

46

46

x = 6234857234503247503427509342579147034917091457093257015340934257093524709435709342570943527045932743590734250975430932457 y = 589258369023486093486094235689081537095123709753120973250973520917090134570349570234957 

40*60*1000/400 

6000

6000

timeit('x*y') 

625 loops, best of 3: 909 ns per loop

625 loops, best of 3: 909 ns per loop

x*y 

367394180509766658895501030664436344529901907226998090136621859173489910\
141757392334581926999855258500352702769438080915735222768105608577713953\
5223872394859098459433166859522089878858692339535983319087299349

3673941805097666588955010306644363445299019072269980901366218591734899101417573923345819269998552585003527027694380809157352227681056085777139535223872394859098459433166859522089878858692339535983319087299349

timeit('x._mul_(y)') 

625 loops, best of 3: 366 ns per loop

625 loops, best of 3: 366 ns per loop

E.conductor().norm() 

79

79

%cython from libc.stdint cimport int64_t def test_64(x): cdef int64_t y = 123507 cdef int64_t z = 509374 cdef int64_t i,j,v j = x w = walltime() for i in range(j): v = y*z return walltime(w) def test_gmp(x): cdef mpz_t y,z,v mpz_init(y) mpz_init(z) mpz_init(v) mpz_set_ui(y,123507ul) mpz_set_ui(z,509374ul) cdef int64_t i,j j = x w = walltime() for i in range(j): mpz_mul_ui(v,y,509374ul) w = walltime(w) mpz_clear(y) mpz_clear(z) mpz_clear(v) return w 

test_gmp(10**9) 

9.3864760398864746

9.3864760398864746

test_64(10**9) 

9.059906005859375e-06

9.059906005859375e-06

4*(1+a)*x+(1+a)**3+(1+a)**3*x+(1+a)**2*x**2 

(3*a + 2)*x^2 + (12*a + 9)*x + 8*a + 5

(3*a + 2)*x^2 + (12*a + 9)*x + 8*a + 5

y = var('y') 

solve(3*y**2+12*y+5-2**31==0,y) 

[y == -sqrt(715827885) - 2, y == sqrt(715827885) - 2]

[y == -sqrt(715827885) - 2, y == sqrt(715827885) - 2]

floor(sqrt(715827885) - 2) 

26752

26752