月: 2024年3月

import numpy as np
import matplotlib.pyplot as plt

def plot(in_curves, func):
    f, (ax1, ax2) = plt.subplots(1, 2)
    ax1.axis('equal')
    ax2.axis('equal')
    for crv in in_curves:
        ax1.plot(crv.real, crv.imag)
        y = func(crv)
        ax2.plot(y.real, y.imag)

in_curves = []
N = 10
for i in range(N):
    crv = np.linspace(-2, 2, 100)  + 1j * (-3*np.pi + i * (3*np.pi-(-3*np.pi)) / N)
    in_curves.append(crv)
for i in range(N):
    crv = -2 + i * (2 - (-2)) / N + 1j * np.linspace(-3*np.pi, 3*np.pi, 100)
    in_curves.append(crv)
print(in_curves)
plt.cla()
plot(in_curves, np.exp)
plt.savefig(ofile)
ofile

(ql:quickload "ironclad")

(defun plain2bytes (message)
  (sb-ext:string-to-octets message  :external-format :utf-8))

(defun bytes2hexstr (bytes)
  (reduce (lambda (x y) (format nil "~x~x" x y)) bytes :initial-value ""))

(let* ((message "hello!"))
  (bytes2hexstr (ironclad:digest-sequence :md5 (plain2bytes message))))

(let* ((message "message digest"))
  (bytes2hexstr (ironclad:digest-sequence :md5 (plain2bytes message))))

import itertools
import functools

@lru_cache(maxsize=1000)
def _calc(i, ind):
    return i * (10**ind)

def Munchhausen(n_digits=4):
    ret = []

    values = [(i, i**i) for i in range(1, 10)]
    for x in itertools.product(*[values]*n_digits):
        a = 0
        b = 0
        for ind, (i, iv) in enumerate(x):
            a += _calc(i, ind)
            b += iv
        if a == b:
            ret.append(a)
    return ret

for i in range(1, 7):
    print(i, Munchhausen(i))

git clone https://gitlab.inria.fr/zimmerma/ecm.git

cd ecm
autoreconf -i
./configure
make -j4
make check

import sympy

N = 100
p1 = sympy.randprime(2**(N+1), 2**(N+2))
p2 = sympy.randprime(2**(N+1), 2**(N+2))

print(f'{p1=}')
print(f'{p2=}')
print(f'{p1*p2=}')

echo 11108239743933603608670116823884965417896735085597594043972249 | ./ecm -c 500 250000

GMP-ECM 7.0.6-dev [configured with GMP 6.3.0, --enable-asm-redc, --enable-assert] [ECM]
Input number is 11108239743933603608670116823884965417896735085597594043972249 (62 digits)
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:4010257250
Step 1 took 77ms
Step 2 took 70ms
Run 2 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2451171661
Step 1 took 70ms
Step 2 took 69ms
Run 3 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2557109986
Step 1 took 70ms
Step 2 took 69ms

...

Run 18 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:504643144
Step 1 took 70ms
Step 2 took 70ms
Run 19 out of 500:
Using B1=250000, B2=128992510, polynomial Dickson(3), sigma=1:2975852467
Step 1 took 70ms
Step 2 took 71ms
********** Factor found in step 2: 2869069113607969685726929882709
Found prime factor of 31 digits: 2869069113607969685726929882709
Prime cofactor 3871722605512470843352197003061 has 31 digits

import numpy as np

p5deg = np.polynomial.polynomial.Polynomial.fromroots([1, 2, 3, 4, 5])
print(p5deg)

def laguerre(f, x0, delta=0.01, max_iter=100):
    ans = []

    f = f.copy()
    while 1 < f.degree():
        n = f.degree()
        df = f.deriv()
        ddf = df.deriv()

        x = x0
        iter = 0
        while not (np.abs(f(x)) < delta):
            numerator = n * f(x)
            denom_plus = df(x) + np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))
            denom_minus = df(x) - np.sqrt((n-1)**2 * df(x)**2 - n*(n-1)*f(x)*ddf(x))

            denominator = denom_plus if np.abs(denom_minus) < np.abs(denom_plus) else denom_minus
            x -= numerator / denominator

            if np.abs(f(x)) < delta:
                break
            if max_iter < iter:
                break
            iter += 1

        ans.append(x)
        f = f // np.polynomial.polynomial.Polynomial.fromroots([x])

    if f.degree() == 1:
        ans.append(-f.coef[0])
    return ans

roots = laguerre(p5deg, x0=0, delta=0.00001)
print(roots)