foxqstm/CheckPrime.py

## CheckPrime.py

# coding: utf-8

# In[29]:


def is_prime(num_in):
    if num_in <= 1:
        return False
    else:
        bound=round(num_in**0.5)
        for num in range(2,bound+1):
            if num_in % num ==0:
                return False
        return True

def is_prime2(num_in):
    if num_in <= 1:
        return False
    if num_in%2 ==0:
        return True
    else:
        bound=round(num_in**0.5)
        bound=(bound-1)//2
        for num in range(1,bound+1):
            if num_in % (2*num+1) ==0:
                return False
        return True

import random

def is_prime3(n):
    """
    return True if num_in is probably prime by Miller-Rabin primary test.
    : param n : int (natural number)
    """
    NumTrial = 100 # Miller-Rabin法での試行回数、素数判定を誤る確率は4^(-NumTrial)

    if n <= 0: return False
    if n == 2: return True
    if n == 1 or n & 1 == 0: return False

    d = (n - 1) >> 1
    while d & 1 == 0:
        d >>= 1

    for k in range(NumTrial):
        a = random.randint(1, n - 1)
        t = d
        y = pow(a, t, n)

        while t != n - 1 and y != 1 and y != n - 1:
            y = (y * y) % n
            t <<= 1

        if y != n - 1 and t & 1 == 0:
            return False

    return True


## ModulesFactorization.py

# coding: utf-8

# In[9]:


import sys

def result_factorization(N,factors):
    """
    return True if non-trivial divisor of N is found
    :param N : int (natural number)
    :param factors : (list of factors of N, len(factors)=2)
    """
    if len(factors) !=2:
        print("len(factors) != 2")
        sys.exit()

    if N != factors[0] * factors[1]:
        print("factors[0]×factors[1] != N @print_factors in ModulesFactorization")
        sys.exit()

    if factors[0] ==1:
        print("自明な約数しか見つかりませんでした。")
        return False
    else:
        print("{0}={1}×{2}".format(N,factors[0],factors[1]))
        return True


# In[1]:


import sys
import math

def is_square(N):
    """
    return True if N is square number
    :param N : int (natural number)
    """
    if N < 0:
        print("N is negative number @is_square in ModulesFactorization.")
        sys.exit()

    sqrt_N=round(math.sqrt(N))
    if N == sqrt_N*sqrt_N:
        return True
    else:
        return False


# In[2]:


import sys
import math
import numpy as np

def GenerateCoprimes(N):
    """
    retruns list of numbers which is relatively prime to N
    :param N :int (Natural number)
    """
    if N <= 0:
        print("N <= 0 @GenerateCoprimes in ModulesFactorization")
        sys.exit()

    Coprimes = list() # Nまでの数と互いに素な数のリスト
    for num in range(1,N+1):
        if math.gcd(N,num) == 1:
            Coprimes.append(num)
    return Coprimes


# In[3]:


def FactorInList(N, List):
    """
    return factor of N in PrimeList if exists
    :param N: int (Natural number)
    :param List: list of int (Candidate of factors)
    """
    a = 1
    b = N
    for num in List:
        if N % num == 0:
            a = num
            b = N//num
            break

    return [a, b]


# In[4]:


import sys

def PositiveInt(N_str):
    """
    return positive integer N for string N_str
    :param N_str: string
    """
    try:
        N = int(N_str)
    except ValueError:
        print("整数を入力してください。")
        sys.exit()

    if N <= 0:
        print("0以下の整数です。1以上の自然数を入力してください。")
        sys.exit()

    return N


## n2a_fixa_PrimeChecker.ipynb

      
Display the source blob

    
Display the rendered blob

    
    Raw
  

              n2a_fixa_PrimeChecker.ipynb
            
          
        Loading

      Sorry, something went wrong. Reload?
      Sorry, we cannot display this file.
      Sorry, this file is invalid so it cannot be displayed.
      
          Viewer requires iframe.
      
    
## n2a_PrimeChecker.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n^2-aが素数になる確率を計算します\n",
      "自然数aの最小値を入力してください:2\n",
      "自然数aの最大値を入力してください:100000\n",
      "各aに対して調べる自然数の個数Numを入力してください:50\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcFOW18PHfmY1hYGBYBkUW2UX2lmFzRcGwKIuggibqlQQTt5uAXjOGuEPCa3L1xqtxTaIxcSVcRFEREsAoggwOgiCrqAygLMqiI8swz/tHd5XdPdV71/RM9/nOhw/T1dVPPdXVU6er6jl1xBiDUkopBZCV6g4opZSqOzQoKKWUsmlQUEopZdOgoJRSyqZBQSmllE2DglJKKZsGBaWUUjYNCkoppWwaFJRSStlyUt2BWLVs2dJ06NAh1d1QSql6ZfXq1fuMMcWR5qt3QaFDhw6UlZWluhtKKVWviMhn0cynp4+UUkrZNCgopZSyaVBQSillq3fXFJRS7jh+/DgVFRUcOXIk1V1RCcjPz6dt27bk5ubG9XoNCkopACoqKigsLKRDhw6ISKq7o+JgjGH//v1UVFTQsWPHuNrIjKDg8cCaNTWn9+sH5eURX15OOQtZyAEOUEQRIxiBB48LHVUqdY4cOaIBoZ4TEVq0aMHevXvjbiMjrinsHdKVqrzsgGlVednsPbNrxNeWU85c5nKAAwAc4ABzmUs5kYOJUvWNBoT6L9FtmBFB4YU7umCyAt+o6mzhhTsiB4WFLOQ4xwOmHec4C1mY1D4qpVRdkBFBYWfrasqu7WUfLVTlZbP62t7sPPlExNdaRwjRTldK1a7HHnuM3r17069fP84++2w2bNgAwLFjx7j22mvp3bs3ffv2ZenSpfZrnn/+eXr37k2fPn0YOXIk+/btq9Hu3XffTZs2bejXrx+9evVi/vz5Cff1zTff5LTTTqNLly7Mnj077Lxz5sxBRAKSddeuXcuQIUPo2bMnvXv3dmVQQEYEhSKK+NcdZ9lHC9XZwj/vOJMiiqJ6bSzTlcoIHg+I1Pznqf1rbVdeeSXr1q1jzZo13HbbbUyfPh2AJ598EoB169axaNEibrnlFqqrq6mqquLnP/85S5YsYe3atfTp04eHH37Yse1p06axZs0aXn75ZaZMmUJ1dXXc/Txx4gQ33ngjb7zxBhs2bOD555+3A1iww4cP89BDDzFo0CB7WlVVFT/60Y947LHHWL9+PUuXLo17hFE4GREURjCCI62bUXZtL6qzYPW1vTl6cjNGMCKq1+YS+MbnkhvVa5VKW0OGQF5e4LS8PDjzzISaHT9+PP3796dnz5488cQTUb2mSZMm9u/ffvutfU59w4YNDBs2DIBWrVpRVFREWVkZxhiMMXz77bcYYzh06BCnnHJK2GWcfvrp5OTksG/fPl599VUGDRqEx+Nh+PDhfPnll1H18/3336dLly506tSJvLw8Jk+ezCuvvOI47x133MFtt91Gfn6+Pe2tt96iT58+9O3bF4AWLVqQnZ3t+PpEZMToI2uk0LI7jnHS+v2suONCJjAhqhFE1jyv8Rrf8i1NaMIoRunoI5XefvEL5xF7lqNHoaoqcFpVlXc039Chzq/p1w/+53/CLvbPf/4zzZs357vvvmPAgAFMnDiRG264gU2bNtWYd/r06Vx99dUAPPLIIzzwwAMcO3aMf/3rXwD07duXV155hcmTJ7Njxw5Wr17Njh07GDhwII8++ii9e/emUaNGdO3alUceeSRsv1auXElWVhbFxcWcffbZrFixAhHhqaee4v777+e///u/WbJkCdOmTavx2oKCApYvX87OnTtp166dPb1t27asXLmyxvzl5eXs2LGDiy++mN///vf29M2bNyMijBgxgr179zJ58mRuu+22sP2OR0YEBfDu3Nu3bk/Zsgu4mhJa0CK219KebWyjF70ooMDFnipVDzRoACedBF98AcZ4Tx2dfHLNo4cYPfTQQ/zf//0fADt27GDLli28+OKLEV934403cuONN/Lcc88xc+ZMnnnmGaZMmcLHH39MSUkJp556KmeeeSY5OTkcP36cRx99lPLycjp16sTNN9/Mb3/7W37961/XaPfBBx/kb3/7G4WFhbz44ouICBUVFUyaNIndu3dz7NgxOx/g/PPPZ02YQGqMqTEteKRQdXU106ZN4+mnn64xb1VVFe+88w6rVq2ioKCAYcOG0b9/f/toKFkyJigAHOQgy1hGV7rGFBQAKqhgLnPpQAcNCir9RfhGD8Du3dCpExw5Avn5sHq1NzDEaenSpSxevJj33nuPgoIChg4dypEjR5g0aVLEIwXL5MmTuf766wHIycnhwQcftJ8788wz6dq1q73j7ty5MwCXX355yIu+06ZN49Zbbw2YdvPNNzN9+nTGjh3L0qVLufvuuwEiHim0bduWHTt22NMrKipqnLY6fPgwH330EUN9R1tffPEFY8eOZf78+bRt25bzzjuPli1bAjB69Gg++OADDQrxmsc8VrISg+FJnqQTndjP/qgT0hrSkFM4hWyyNZlNKYDWreHaa+Hxx73/JxAQAA4ePEizZs0oKChg48aNrFixAiDikcKWLVvo2tU7vHzBggX275WVlRhjaNSoEYsWLSInJ4cePXqwa9cuNmzYwN69eykuLmbRokWcfvrpAPYF55tuuilsP9u0aQPAM888Y0+PdKQwYMAAtmzZwvbt22nTpg0vvPACzz33XMA8TZs2DRgJNXToUH7/+99TUlJC586duf/++6msrCQvL49ly5Y5BqFEZURQmMc8VrDCfmwwbGOb/dhKSAMi7tw3sIG3eMvOXYjltUqlnTvugPXrvf8naOTIkTz22GP06dOH0047jcGDB0f1uocffpjFixeTm5tLs2bN7B31nj17GDFiBFlZWbRp04Znn30WgFNOOYW77rqLc889l9zcXE499VT7dM3GjRs566yzwi7v7rvv5rLLLqNNmzYMHjyY7du3R9XPnJwcHn74YUaMGMGJEyeYMmUKPXv2BODOO++kpKSEsWPHhnx9s2bNmD59OgMGDEBEGD16NBdddFFUy46FOJ3nqstKSkpMrEV2bud2DJHXs4giSil1fG4f+9jCFpaylIMcjOm1StUHH3/8sf2NOVNdfPHFzJ07l7wEr42kmtO2FJHVxpiSSK/NiCOFaAIChE9I28UuXsF5+Fik1yql6ofXXnst1V1IuYzIUxCiuxdIuIQ065pCE5o4Pq/JbEqpdJARQWEQgyLOEykhTRByyOEcztFkNpW26tvpZFVTotswI4LCeMYzmO8vWgnCqZxqPy6iKGIyW1Oa0pvenMEZjGZ0TK9Vqj7Iz89n//79GhjqMauegn8mdKwy4poCeAODNQLpEi6hN725h3sYxjAu5MKIr9/NbhawgG50ow99eIVXGMtYziSxtH6l6oq2bdtSUVGR0L34VepZldfilTFBwb/+wVzm8iqvUkAB/+SfrGZ1QK6BUx5CYxpTRBF/4k8c4hCNaRzVMhPJZ0hmPkRdzK1w6hNQ5/qZKXJzc+Ou1qXSR0YMSS2nnJd5mWpC3+Ewl1wmMAHwBg3/Ggq55NKFLmxkY8BIpmyyuZRLHXdaVnGe4HaiPdWU6OvdaitZnPqURRaCcILvb2me6n4qlS6iHZKaEdcUFrIwbECA7wvnhCqqExwQAE5wImSxnUSL8ySzuE9dLBTk1KdqqgMCAqS+n0plmow4fRRtDkG4+ULlOsRahCfRvsSTD1EXCwXFsmzNAVGq9mTEkUK0OQRFvh8noXIdYi3CE0tfEnm9W20lSyzL1hwQpWpPRgSFEYwgK8KqWrkGIxhBTtABVC65dKd7jcCQTXbI/IREi/Mks7hPXSwU5NSnLLLIJrBoSKr7qVSmyYjTR9ZFyhf5/m6LueTa57SDR7kc5Siv8AoGYz93CqeQRRbb2U4llQD0oEfIC6DW9Nd5ncMcjnkkjTXfHOZwghMJjcQJXv+6MKonVJ/g+4JGdaGfSmWajAgK4N0J+QeFGczgbu5mCEMYx7iAeQcykD70Idf3A7CWtaxnPdOYRmMacx/30YEOEZeZyA7N6nNzmnMbiVVY8uBhLWvpSU9KiDgAoVZ48FBIIUUU0ZKWAdOVUqnh6ukjERkpIptEZKuIhLyFqIhcKiJGRGplb3U5l5NPPv3pz8d8XOP5Siq5l3sp4/uhrwUUcAqn8Cqv8h7vcRM30Yc+EZdVRhnP8EzE+UIZzGBGMSru1/vbzGb2sS/yjLXoaZ5mFasCpr3LuwEBXClVe1w7UhCRbOAR4EKgAlglIvONMRuC5isE/hOoWazUJS/5fizzmMdGNtoJU+dyLgBVfF+DNoccmtCEjWxkK1upoIJiilnHOg5wgIY0RBD71FK0IiWVjWd8gmv7vSu50v5G7r9cQQJOlVnL9y9MJAiDGBR1f8Ktl/9zOeTwDd8EvPZVXgVgEpPiaj/ZyilnPvP5ju8A7xeEMYxJaHmx9r+21jfZ/UpWv+tC8mUq+pCKZbp5+mggsNUY8wmAiLwAjAM2BM13H3A/cCsp4l+A5wAHeIM3aESjgKzlfPJpTWvO4Awa0IC/8Bc2s9nOf7B2GE6yyKKc8hobMziBy6lgz8M8TD/6cTZnJ7yez/Ec53AOJ3FSwHKt4bb+y/+Mz2oUJrIeRwoM4dYLApMDq6iinHK60MVe5ylM4Su+iqv9ZP/BOCU+VlLJHObEvbxY+19b65vsfiWr37W5vetSH1K13m6ePmoD7PB7XOGbZhMRD9DOGFOnbmJ+nONkk01PetrT9rKXZSzjJE6iHe0AIibEWaqpdkzAiiaprIIKlrM8ntWooRvdaElLx+UGL39liAO3UNP9hVuvUElr/uvcjW4BNzCMpf1kC5X4GC5xMZo2Y+l/ba1vsvuVrH7XheTLVPQhVevtZlBwGthvZ4CJSBbwIHBLxIZErhORMhEpq62bdR3iEKtZbT/OJ59WtGIOc3ibt2NuzykBK5qksrM4i4tITsk965pCpGSwAxwImawXTcGicOsVzTovYxl/5a9xtZ9s4dqMd3nJSnhM9vomu1/J6nddSL5MRR9Std5uBoUK8H2l9moL7PJ7XAj0ApaKyKfAYGC+08VmY8wTxpgSY0xJcXGxi10O5H9NIY88WtCCHexgKUtr5DJE4pSAFU1S2RjGBByxJOJKrqQf/SImgxVRFDJZL5qCReHWK5p1foM32FDjLGN07SdbuDbjXV6yEh6Tvb7J7ley+l0Xki9T0YdUrbebQWEV0FVEOopIHjAZmG89aYw5aIxpaYzpYIzpAKwAxhpjYrvbnQtyyaWQwoAqa7nk0oIW/Af/wVSmUkVVxIQ4SxZZjglY0SSV/YE/8G/+HeeaBHqe5ymn3HG5wcsPVZgomoJF4dbLKTkw+P2ZylQmMjGu9pMtVOJjuMTFaNp0SpB0KxEyln7FspxI8yer33Uh+TIVfUjVert2odkYUyUiNwELgWzgz8aY9SJyL1BmjJkfvoXakUUWAxnIWtZSSSVFFHEhF9KZzjSkoT3fV3zFe7xHf/pTSCEAfenLJjZRSWXI0UfhRqpEk1S2m90sZznncE7C69qNbhRTbLe/gAUBI3/8l2/NE8/oo+D2G9OYi7jInn6EIyxkIUc4QiMacR7nBaxzZ99PpPZrI8nNatP/gl+io488eDjGMd7gDY5wJGL/ren/4B9UUeXa+gYvpylNGcnIiP2yLsQH98v6/1Vetf+24um3NX8qRx+lIgE0VUmnriavGWNeB14PmnZniHmHutmXUJrQhPGMpx3taEELOtCBb/iGmcxkHOMYwhDAe/qoKU15jufoSldu4RYa0Yge9KCSyhrfoLeylV3sIoss3uXdsH9YL/IiHejAz/hZjeeHMpT2tE/Kum5iE61oZS+3JS05xCE+4AO60a3GOoxnPBdwATnkUEBBTMvy4KEpTXmap/kRPwpI9OtFL6qoogc9eIAHagTSxSzmMz7jx/w4bPsFFLCLXZzP+TH1LVYePGSTzbu8y1VcFVUtjUh605tKKjmd0zmZk6PqwzrW0YUurhZ28uBhDWvYxCZu4RbyyIs4/yIW0YIWjtvLg8cezu1f7TCefqU6qdE/AbaUkGlXSV/mPOZxlKO1tsyMyWgOxbpw+jIv05e+dKCDPc3/mkJDGtKGNnzIh3zJl+xkJ2dzNmWUsZe9DGIQt3N7VBdig/2aX4c8nTOSkXGslbMruZIWtLAfL2IRm9lMJzqFHI30G37DSZzENKbFvLwmNGEgA+0jK8vXfM0CFlBEER3pWOMc6WIWR9X+RjayilWuBwXLZ3zGTGYym9kJt/Ut37KQhRRSGFVQAO/nMdoRb4m4lmtjmn8YwwKOqoP9nb9zGqfF3G4wa6eYjPc/XqlY9j3cU6vLy4igYCWAODnIQfvD9qHvx7LA9+OvAQ04ylG2s53GNA64iVsOOTV2rmdwBlvYQimldpKYf7JYd7qzhjUc4UiNvlnz5ZPPOMbZY7+dDqP9pwe30YlO7Gc/BzhAAxrgwcOHfGjnVnzKpxRS6JjQBt7TPcEJXME605mpTA2YZjB8wRc8zuMc4pDdpvVNWxAqqLCH+Fqu4zoOcjBgmtN672MfVVRRSmnA6btQ70sih9/b2W7/PpvZAe07JbZ9xmcBp978t0ERRfa35jnMYTGLIyZ9dac7n/AJm9nMAhYEfH78Ey+DkwTjTbqL5X1zqoPh3w54j1L937dYlm9t2+A2o+mfUzuhTmVFm3AZ6+comiTQupCcZ0n7ymtOFb6SJZtsBjOY1rSmhBJXl5VLLv3pz2pW16ig5jQ9HtlkO/5xZ5GF8f2EExwYFrKQJSwJuaxRjKIlLWlOc/u0lpNQVdogdK5IuPcr1kpu5ZTbNyYMbv993k/Kt3f/fiXyOfKvIOhUbTBctUDLAzzAHvaE7F+w2cymOc25jusCpsdb8S/S+kdboS9SO5He80jVGKP5HM1jXkASqGUwg+3AEOl9+jW/poqqhI9StPKaT7hErUSd4ATllNvtu7ms4xxnJSsdk1mcpscj1Le9aqqjOi22jW0Bj4PvaRS8rLd5m2d5lg/4IOC5N3mTP/JH+3GohLdwO+Nw71c8yVNOFeFWsjJpp3MiJX3F2k4iSXdO2eTh3jf/G0f6izf5KtL6R1uhL1I70SbaJZJEFk0SaKT2G9Ag4nKSKe1PH7md6FFJJctYxhCGuL6sRBLKUiH4nkbBDnGILnSpcU1hKUsDHsf7vsZaLS+UUPMn+32PlPQVazvxPu9/LS2a1w1jGPnkRz1/ov2L9nXRtJNIol007UfzNxtp2XdwR8TlJFPaHym4neiRS64dyaNJ7EpEIgllqRB8gTlYIxrxBV9wmMMB06/neq7kSvtxvNsw1mp5ocRajS9ekZK+YmknkaS7WJOmjvp+Em0n2uejfV007USTaJdIElk0f7N1ITnPX9oHhXCJWomyrikMZSjgvSW3W8vKJZdBDHJMehrEoKQsN7jqmcU6hxtJcG7BAAaEXdb5nM94xtOb3gHPncqpAbclD1WlLVzyYKj3Jd7kKaeKcIMYFHUCYySRkr5ibSeRpDunbP1w79tiFjvenyvWJD3/14Vb/2gr9DktP9RrIiVcxvs5iiYJNFL7pb6f2pL2QcGDhwlMiPitFWL75teEJoxjHKdwin2axFqWFeHzyQ8oHhOsiKKwO05LPvlMYALjGc9ZnBXwemv6BCY4jikXhFM4JaCtwQwO+BBmkUU/+nEpl9p9t96LQgoZwxgu5/KwfXQafdSJTgA1his2oQndfT/P8zxrWBPw/Gu8xh/4g/3Yg4fxjLe3YRFFXMZldg6JtQxr/YPfl+D3K57kqQEMsHey/u1PZGLAqZMCCpjEpIAb+gkSEDCLKKI//e0jzOB+WZ8ja3lNaFKjPet1gxkc8G3XaseDh8u4LOC9L6Ag4kVmgNa0Dngc6X0r8P04vW+XcIndh2jf/+C/I/91aEpTLuOygM9qqHY9eBjHOHu0W0Ma2v0M9Z47rXNwf4CoP0fjGV9j2/lfZPZftv/29m+/KU0jLieZ0v6aAnxf4espnnJ8/nzOZwlLEIR2tONrvg44pZFNNudwDgc4YO/ACikkn3zKKONLvuQczqkRzY/4fgB+zI+ppJJVrCKHHAYwgJ705AhHWMUqTuIkDnGI4Qy36wlcxVU8y7NczdX2DrY73VnGMn7CT+hCl4B1zCGHRSwim2x2sxuAX/ALCihgLnNpTGMGM5g2tOFTPuULvgCgLW05lVPx4KGYYj7kQ4oppjOd+R2/YznLuYVbyCKL53meRjTil/ySO/HmIU5nuuPooWY0YxjD+JRP2cc++tKXaqrpRS8e4zE7LyT4MPkd3qnRVh/6MIc5jGIU53EeAJ/zOQCXcAmDGMQDPMBJnMQP+WHA+3KYw+xnP5dwieP2j0ZrWtsXbn/JL+0dcy962cOWm9CEX/Ere7nWqJPf8lsAPuAD1rOeq7iKvexlNavpRjemMKXG8jx4yCOPMsq4givII89u73Iu5xu+sbPcN7KR93iPSUwK2Dl78LCHPSxhCcN9P9Gw3r85zKGa6ohfCIYzPOTF0JM4ie/4jmyyY/q268HDUY5SQAF96MNi389BDto7y7a0ZTe76UnPkEe5HenIN3zDJCbhwcMWtvAyL3MN19QIfh489u3wr+CKgOd60pONbGQ84+lO96jXA7yBwdp2pZQ67uSdtrfldm6PaXmJSvshqeB8T/zaYh2+VlEVMPY/lLGMpRvdeIVX2MrWgPmLKKI5zdnOdlcuLueQE3CRsR3t2MEOe6hqEUX0pS+HOMQnfGLnEoxmNOdybthciSyy7Dba0paP+Mh+rgc9uIqr7Pl3sIM5zOFLvrSndaITX/IlxzgW1cic4JwQJ/7j9kP1PZtssshyXGYBBfSgR0CFPusIzJo/hxwmMrHGMvz7Ze2w/DnNm002PejBOtYB3ydSzWc+y1nOrdxa48h0LnN5n/cd1z84zyH4PXMa1w818wOOcIR88h3H9Ft5Pf7bxKkt/+da0CJgJFseeVzCJfbtMqz1nsUsDnOYJjThEIdoSEOqqAq4HUl3uvMhH3KCEzSkIdVUc5SjNKAB2WQ75rU45XZ0pSszmUkPetCEJo55B+FyiKys5GhyRdzKWYh2SGraB4VyyutVaccsslISvIIJQne6s5WtEXfCOeRQQkncuRL+O06AJ3myxvBWt2STzQAGJCXPI9IyVrHKcdhv8Jj3SOPre9ObRjQKGOe+kIXcxE01bsHxFE+xla1JWY9Q+QFZZHEyJ3M910edY+HUVjiCcA7nMJrR9rTXeI13eTcpX5CsvBOnnIJsshnHOJaz3D66DtaZznzO51HlEIXKcdjFLt7hHdaxzvH90zyFJKnNQhzJUBcCAniHzG1kY1Q7yiqqEsqVqKIqYDvVVkAA77j9ZOV5RFpGqB1g8Jj3SOPrd7CjxjnpUkod78m0l+TVHwmVH3Cc4/ayo82xCJcF7cRgeJu3AzL/P+KjpB0xW3knTk5wgn/xL37BL0K+fhvbos4hCpXj8Df+xgd84OpnMRppHxRqsxBHuonlDy7RP85UbqfayPOItAz/9Y9mHL9/LsfHfMxTPMW3fFtj3uDbhbihmmr7Jn1ub8e7udv+PdnLCreNDnCAZ3k2aW069X0sY2Nu3w1pHxRSNdY30yQ6Zj+V26k28jwiLcN//aN5L97kTfv3LWxhK1sd70sV6S6nydCQhnzN14C72zH4PUz2dgvXXlOasp71SWvT6X3qTvc6kbOQ9kGhNgtxJEOyxr0nShBO5/SoxsvnkJNQrkQOOQHbKVwthWTLJjtpeR6RlhFqhEzwmPdI4/T70CdgaPKpnEpLWjpmFSdyu+pgofIDqqm2b2kSbY6FU1vhCMIQhgScVz+Xc5MWGKy8EyfZZDOMYbSlbcjXd6azY65BLLkyO9lJBzo4vn+1eSRdN/ZALvLgYRKTauXbYD75CX2baUhD+wPhdH64iKIaw+iSyRopVUQRl3M513CNPfzT3+mcHvB4vO9nAhMchyYGB7rgC2xd6RowzelbUSc6MZKRjju+cMK9/9a4/XB5HuF2XAUUMIhBAX2yLrz6s5bhdA69MY1rXHQMHjMf7EquZAxj7Md96cut3Or4mbmGa7iUSx3bsfIfovkiYuWGOOUHNKOZvWyr79bQWP/PtLUtGtPYbiv4/bVup+4vjzwMhuUsDzhFNsr3Y31uG9IwIFmtgIKAPCD/551yFpwKSFmfkYEM5CZuqpEMZ+UdTGUqE5hg51RYuQax5Mo8y7OsYQ0jGRmQE+NfAbI2ZEyeQh55cZ0TjKQpTamkkh/zYx7jsRrPW+cUr+RKvuVbNrOZnezkEIdqzHuMY4xhDPOYF3DfoJ/yU/sPZTvbeZzHAe997//CX+hABz7lU8A7jLSKKjtPAaCEEsoooz3t7bH9wYYxjAu5kM/4jFWsYj/7+ZIv7doG7WnP+ZzPMzxTow2rCJAHD1/xFatZbd9UbSQj2er76UEPTuZkutHNvv1xF7pwGqcFtLea1YA3f6Qf/TiJkzjGMe7kTkYxitWspgtd2M3ugFtaW67mav7KX7mYi+lGN57hGRrRqEa/7+AOe0flwcMCFnCMYxRRRDOa8VN+CtQsqFJAgZ2jcZSj9gXKHHK4j/v4nM8xGB7lUbtt8I51X896+tLXvkV7K1o57iA8eGhAA1awgq/4in3ss597n/c5zGGGMQyADWzgbd7mKq6iEY0C2nmLt3ibt+3HU5lKG9rwFV/RkpYB+Q8jGUlzmrOTnSxjGa1oxR72cDu3B4ytD+5vAQUBAdXKd9nOdtrS1v68Dmc4i1jEMIbZbRzjGIUU8h3f0ZnOFFFEJZXcy712e/dyL8tZznzmcx/32UcLO9nJMY5xCZfYr7W21xCGMI5x7GOffRRzF3fxv/wvO9nJT/gJbWhT430fyMCAIbzWdj7GMZ7jOSYzmV70Ygtb+IiPuJiL7aDnwcMqVvEJn3AN19jte/DwEi9hMPyUn9KMZjWWCzCRiXzERwxkIF/zNY1pbN8toTal/ZGC5TM+c6XdSio5znHHgABwK7cC8BzP8QqvcDInM4pRFFPMndzJDGYEzH86p3Mrt/IjfsQ7BbUIAAAa00lEQVQoRgHwOI/bO6aOdKSQQk7jNDuRzj+zdwc7atzQrIIKhjGMz/mcYood+/kRH1FKKY/yKGWUsZjFzGWu/fwkJnE6pzOFKVzABcxkpl0cxv/bXmMa05KW9h/1m7xJc5oD0IY2/IAf2N+2ruZqDnKwRi2Jn/NzRjOaJSzhQR4EvN/IWtCCd3mXPexhC1vsgHAbtwX0YQELKKGEtrTlC75gP/vtZU5gAqWUMpvZNY4ibuZmbuAGDnCA7Wy3L2rezM0B81VSyV/4i90vi7Web/GWHRCsHTfAGMZg1b/uSU8APuETx+0B3gu4m9nMOMYxjWm0ohXDGMZc5rKIRfZ829jGp3zqeE0h+L19kif5hE94iIfskUmXcRk/5IesYhVb2EJf+jKc4fbnL/g22utYRyml9iixL/iixiin7WxnAQvsb/b55HMe5/ErfhXw7f11Xucd3uFlXmYnOwHvt/jgI+IzOTOgQBTAX/kri1nMy7xMBRUAPM3TAHalP2uHfREXAdh5HCc4wau8SimlAUdwE5hgFwPyT3Y8znE2spEyyqimmj/xJ8e75P6EnzCLWQF3EQDsI26rn0660pVLuIRccu3aG6mQ9nkKEPqe5m6IJkEtHVjj7mMliGMxokia0jTkSBqnOhDhakNYf8hW8Ztd7ApZPMjantbrTuM0NrEJ8B7aD2c4c5hjz19EUY0/5tnMrjF+P598jnCEC7iAYoprJCt9xmcBn9lQ6+OvP/3ZwIaQ65KoZHy288jjGMdqJJqFEmq988ijM53ZzGZOcCKqvhVRRCGF7GBHwHSrFjVgf6OH6D4fwUcW2WTTgAZUUkkuuVRRFbJfodbNWhcr0c0/zyqRZDZNXvOpzYCg3GV9K69LQTeXXEYxioUs5BjHHPs2iUkhM+p70IMtbEn52PRMF2vSaG19+XPqVzyFokCT12yhElJU/RNN9bfadpzjzGc+Rzkasm+hCt6A93qABoTUizVptLY+h079iqdQVCzSPijUtZ2IyjyaQKmSzc3PVNoHhbpagEZlDk2gVMnm5mcq7YNCqIQUVf+I76cuySWXcYyjgIKQfRvBiJDP9aCHq4lzKjqxJo3W1ufQqV/xFIqKbZlpLrjIhapd1jBNiyAxZbJamtIUg6lRpc1JtO0LQnvaO/7hWX/0Tn/8/kVuxjCGV3iFSirtESMQWGzdg8ce3hnsG75hAhPsPIAiipjEJAYxKOYdT3/616kAE67/DWmYUPZ+Hnn0oEfY7RSsiCLH5EcrMc8pydUpu76AAs7mbPrRj4EMDLm8SNsi0ue0AQ3sfgUXkErGrbRDSfugADhmKtaWEgIv9ueRZ9cf7kSnGkW5/4v/sn+va9+Kz+VcpjCFMYxhJjPDznsGZwDYSWzgTV76Lb/lfM4HvDvL4FsHON1KYCITuYVbaEELx9tAB39rKqSQMzgjoGKXlUjkXwxoBjMYz/iAi3lWkqB1Lern/JzJTA5o32BoT3tKKaUf/ezprWhlz2vVELAK21gJfqdzOjdwA/3pD3gLBVmV5c7jPEopxYOHi7jI7sNUpkZ1s7ShDLXH43ehC61p7bjjcTr1cDEXB/zel74ByV2zmR0wT7AxjLE/69Y6B1/Pm8UswDuc+S7uogc9QrbnlJ1tmcQk7uVedrHLXsZoRkcscF9Kqb1M/8/Zf/FfdoW1u7jLnn45lzOVqfa2s/ySX/Jv/s0a1nARFzGGMRRQwPVcb+coTGAC93EfeeTZORP+ZjOb27mdPvQJ+DufzGS7qFFf+tr9yiGHIQyxPx9uSvshqaDDUtX36kq9ikzUmc6UUGLnZESTe9Ga1pzLuQGFb+JhFer5mI9Zy9q420k1az00T8GP5ikolVmiCR7RqIt5LvHIIovLuEzzFOKleQpK1W/JCAhQN/Nc4lFNteYpJCIdPgRKKeVP8xQSUNcu1iqloldX6ovUNZqnkADNU1Cq/qqmOq4hzE7qYp5LPLLI0jyFRGieglLfE8SxjkBtCM5ZiTanIviaQrjhquHavJzLKaQwqmXWVbnkxnWRORZpHxQgtXkK4dS32x+ES9RJlngqywWPA88hh050cpzXSi4LxT/voLbMZjY/5IeAt0hSMKvYTzRCFXCxGExA7ojFv85GcC0AS1/6hmx3JCM5kzPDLvs2bqMRjcgii4u5mC50CTmvf5JZcLWzH/CDGjv/8zmfX/GrsDcX9OChO90dn+tGN7LIsnNoLP61HyyhEhGtWg2WFrSgOc1rVCoMFipQdac7AxkYEAAmMtH1PAVXg4KIjBSRTSKyVURKHZ7/mYisE5E1IvKOiITOZklD9e1Gaf73jffnX+SnG90SWoZ/xbhoWVXnLFVUORavKaAgZJEh8Fayi1TT2H9dnTi9vglNHMuaWkop5e/8HaDGvf4Bu3JZNELt0P05jff3L5Kzi12Or7OqxTl5kzeppDLisu/gDtrTng1sCLstrOJArWjFTGYykpF2MqJ/XQrLEpbwG34Tsj2rJGbw0YplM5upppolLAmYbtUM8c9Qf4M3HNvYx76AL3q96MUYxgRUWfSvrW25juvs363CVQAb2cj7vG9XKbyQCyMGmGRwLU9BRLKBzcCFQAWwCrjCGLPBb54mxphDvt/HAjcYY0aGazee5LVyygMKVSilVH1kFd5xM3nNzSOFgcBWY8wnxphjwAvAOP8ZrIDg0wiSP35UA4JSKl1UUskc5thHD25wMyi0gYBj4QrftAAicqOIbAPuB/4z2Z1wM8lDKaVq2wlO1NvkNaexXzWOBIwxjxhjOgO/BH7t2JDIdSJSJiJle/fudZolpPp23l4ppSKpr8lrFRAwlKIthLiC5fUCOA8TMsY8YYwpMcaUFBeHvjjlpL6N8FFKqUjqa/LaKqCriHQUkTxgMjDffwYR6er38CJgS7I74WaSh1JK1bZssl3dr+VEniU+xpgqEbkJWAhkA382xqwXkXuBMmPMfOAmERkOHAe+Bq5Jdj+sq/R6sVkpVd8lcuvsaKX9rbMtpdRIk1BKpUgxxRzgQNhkMxVoOMM5j/Pirq5XF4ak1iluZwEqVRfV1ds67GVvnQgIThnksfCvlBjsbu5OqO1gX/N1rZRbzYgjBc1VUEqli850ZipTY36dHin4aEBQSqWTbWzjSZ50rf20DwqavKaUSjfb2OZa22kfFDR5TSmlopf2QUGT15RSKnppHxQ0eU0plW4609m1ttM+KHjwMIlJqe6GUkolRbyjj6LlWkZzXeLBwxa28AEfpLorSikVt9nMdn0ZGREUHuAB9rAn1d1QSqmElFJKK1oxnemuLSPtTx9pQFBKpZM97OEBHnCt/bQPChoQlFLpxs39WtoHBaWUUtHToKCUUsqW9kGhFa1S3QWllEoqN/draR8UpjNdA4NSKm24PfooI4akTmc6L/GS5ikopeq13/Abslz+Lh9z6yLSSESy3eiMm/LJT3UXlFIqITvY4foyIh4piEgWMBn4ITAAOAo0EJG9wOvAE8aYLa72MkGaq6CUSgeP8ig55DCRia5Vk4zmSGEJ0Bm4HTjZGNPOGNMKOAdYAcwWkR+50rsk0ICglEonVVTxEi9RTrkr7UdzTWG4MaZGMVVjzFfAP4B/iIj7hUPjpAFBKZVuDIaFLHTlaCFiUAgOCCLSCr4/QW+M+dwpaCillHKPWwXEor7QLCJjRWQLsB1YBnwKvOFKr5RSSoXlVgGxWEYf3QcMBjYbYzoCw4B3XelVEmmOglIq3QjiWgGxWILCcWPMfiBLRLKMMUuAfq70Kok0eU0plU5yyOFyLndt9FEsyWsHRKQx8DbwdxHZA1S50qskm8505jOf5SxPdVeUUipuk5jkWjCwiDEmuhlFGgHf4T26+CHQFPi77+ih1pSUlJiysrKYXjODGZzghEs9Ukqp2hVPSU4RWW2MKYk0X8TTRyIiAMaYb40x1caYKmPMM8aYh6yAYM1TF2lAUEqlm21s40medKXtqJLXRORmEWnvP1FE8kTkAhF5BrjGld4lgQYEpVQ62sY2V9qN5prCSGAK8LyIdAQOAA3xBpS3gAeNMWtc6Z1SSqlaFU3y2hHgj8AffZnLLYHvjDHuZE4opZRKmWhuiJcP/AzoAqwF/myMqRejjgCyydZTSEqptNOZzq60G801hWeAEmAdMBr4b1d64pJZzCKbenenb6WUCime0UfRiuaaQg9jTG8AEfkT8L4rPXHRLGZpnoJSKi3MZrar7UdzpGDf7C7W00YiMlJENonIVhEpdXh+uohsEJG1IvJPETk1lvZjoQFBKaUiiyYo9BWRQ75/h4E+1u8icijUi3zV2R4BRgE9gCtEpEfQbOVAiTGmDzAHuD++1QjPrfuOK6VUbSul1NV9WjSjj+I9IT8Q2GqM+QRARF4AxgEb/Npe4jf/CiDpxXrKKedFXkx2s0oplTLWPs2NW164WQG6DQQUFK3wTQvlx7hwK+6FLEx2k0oplXJu7dtiuSFerJxufeF4oyVfOc8S4LwQz18HXAfQvn17p1lCcqsQhVJKpVLKi+zEoQJo5/e4LbAreCYRGQ7MAMYaY446NWSMecIYU2KMKSkuLo6pE24VolBKqVSqC0V2YrUK6CoiHUUkD5gMzPefQUQ8wON4A4IrxZTdKkShlFKpVBeK7MTEN3z1JmAh8DHwkjFmvYjcKyJjfbP9DmgMvCwia0Rkfojm4ubBwyQmJbtZpZRKGTfrKrh5TQFjzOvA60HT7vT7fbiby7d48OgIJKVUWqgLyWtpoQUtUt0FpZRKyEhGur4MV48U6oq7uIujOF7DVkqpeuNN34+b9z5K+yMFDQhKqXST6spr9ZoGBKVUOnKr8lraBwWllFLR06CglFLKlvZBoQENUt0FpZRKulRWXqvX7uEeDQxKqbSS6spr9d493MMf+AO72Z3qriilVNzGMIazOMvVZaT9kQJ4y3FqQFBK1Xev8iqllDKPea4tI+2DwixmcZjDqe6GUkolzQpWuBYY0j4oaEBQSqWjlax0pd20DwpKKZWOjHPNsoRpUFBKqXpIHItbJi7tg0IhhanuglJKJd0gBrnSbtoHhRnM0MCglEorgxnMeMa70nZG5CnMYAaP8zjb2Z7qriilVNwmMpEBDHB1GWl/pGBpQpNUd0EppRKyl72uLyNjjhROcCLV3VBKqYS87fvJJptZzHJlGWl/pKABQSmVbk5wghnMcKXttA8KGhCUUunIrX1b2gcFpZRS0dOgoJRSypb2QSGb7FR3QSmlks6tfVvaB4VZzNLAoJRKK26OPsqIIamzmMVzPMda1qa6K0opFbchDGEc41xdRsYEBb2FtlKqvnvP91NIoQ5JjZcGBKVUujnMYU1ei5cGBKVUOnJr35b2QUEppVT0NCgopZSypX1Q0FoKSql05Na+Le2DghbZUUqlGzdHH2XEkNQZzGABC/g3/051V5RSKm4DGMBEJrq6DFePFERkpIhsEpGtIlLq8Py5IvKBiFSJyKVu9qWSSjebV0op161ilevLcO1IQUSygUeAC4EKYJWIzDfGbPCb7XPgP4Bb3eoHaE0FpVT6KKW03hbZGQhsNcZ8Yow5BrwAgfnZxphPjTFrgWq3OqEBQSmVbuprkZ02wA6/xxW+abVKA4JSKh3VxyI74jDNxNWQyHUiUiYiZXv3ul+4WimlMpWbQaECaOf3uC2wK56GjDFPGGNKjDElxcXFSemcUkqpmtwMCquAriLSUUTygMnAfBeX50hrKSil0lG9K7JjjKkCbgIWAh8DLxlj1ovIvSIyFkBEBohIBXAZ8LiIrE92P7TIjlIq3bg5+kiMies0f8qUlJSYsrKymF/3PM/zIR+60COllKodOeQwk5lxvVZEVhtjSiLNl/a3ubBoQFBK1XcmvrE6McmIoODWYZZSStWmE5yglFJX92lpHxS08ppSKt1o5bUEaEBQSqUjrbymlFLKdRoUlFJK2dI+KGiBHaVUOtLKa3HSymtKqXSjldcSNIMZ/JW/soENkWdWSqk6qgENuId7XF1G2h8pAMxjngYEpVS9d5SjlFLKAzzg2jLSPijMYx4rWJHqbiilVNLsYY9rgSHtg8JKVqa6C0oplXR72ONKu2kfFGrjXiFKKZUu0j4oiGMBOKWUUk7SPigMYlCqu6CUUknXilautJv2QWE84xnM4FR3QymlkqYVrZjOdFfazog8hfGM5xzO4Xf8LtVdUUqpuPWkJ1dxlavLSPsjBcu3fJvqLiilVELcOmXkLyOOFAD+yB9T3QWllErIEt+PmzWaM+JIoZTSVHdBKaWS5gQnXLv3UUYEBaWUSjcnOOFKuxoUlFJK2TQoKKWUsmlQUEqpeiibbFfazYigMJvZqe6CUkoljZujjzJmSOqN3MgjPJLqbiilVNyGMpSRjHR1GRkTFF7jtVR3QSmlErLU9+NmBbaMOH30JE/yGZ+luhtKKZUURznKXdzlStsZERS2sS3VXVBKqaQ6ylFX2s2IoKCUUio6GhSUUkrZMiIodKZzqruglFJJ1YAGrrSbEUFhKlM1MCil0oabo48yZkjqVKYGPP6QDznOcUooCbiLag45zGRmvb+zqgcP5ZTH/foWtGA/+5PYI6VUos7jPEYxytVluHqkICIjRWSTiGwVkRp7WRFpICIv+p5fKSId3OyPv2/4Jq0L7+xhT0KvT+f3Rqn6Ko8815chxhh3GhbJBjYDFwIVwCrgCmPMBr95bgD6GGN+JiKTgUuMMZPCtVtSUmLKysoS6ls55fyDf1BFFbnkcpzjgX1HMLjzviilVCJyyGEiE/Hgiel1IrLaGFMSaT43jxQGAluNMZ8YY44BLwDjguYZBzzj+30OMExExMU+UU45c5lLFVUANQICoAFBKVVnVVHFS7yU0OnhcNwMCm2AHX6PK3zTHOcxxlQBB4EWLvaJhSx0DARKKVVfGAwLWehK224GBadv/MFfwaOZBxG5TkTKRKRs7969CXXqAAcSer1SStUFbu3L3AwKFUA7v8dtgV2h5hGRHKAp8FVwQ8aYJ4wxJcaYkuLi4oQ6VURRQq9XSqm6wK19mZtBYRXQVUQ6ikgeMBmYHzTPfOAa3++XAv8ybl359hnBCHLJdXMRSinlKkEYwQhX2nYtT8EYUyUiNwELgWzgz8aY9SJyL1BmjJkP/Al4VkS24j1CmOxWfyzWFfuFLOQAB3T0kVKqXol39FG0XBuS6pZkDElVSqlMUxeGpCqllKpnNCgopZSyaVBQSill06CglFLKpkFBKaWUrd6NPhKRvcBncb68JbAvid2pD3SdM4Ouc2ZIZJ1PNcZEzP6td0EhESJSFs2QrHSi65wZdJ0zQ22ss54+UkopZdOgoJRSypZpQeGJVHcgBXSdM4Ouc2ZwfZ0z6pqCUkqp8DLtSEEppVQYGRMURGSkiGwSka0iUprq/sRCRNqJyBIR+VhE1ovIz33Tm4vIIhHZ4vu/mW+6iMhDvnVdKyJn+LV1jW/+LSJyjd/0/iKyzveah9wuixotEckWkXIRec33uKOIrPT1/0XfbdkRkQa+x1t9z3fwa+N23/RNIjLCb3qd+0yISJGIzBGRjb7tPSTdt7OITPN9rj8SkedFJD/dtrOI/FlE9ojIR37TXN+uoZYRljEm7f/hvXX3NqATkAd8CPRIdb9i6H9r4Azf74XAZqAHcD9Q6pteCvw/3++jgTfwVrYbDKz0TW8OfOL7v5nv92a+594Hhvhe8wYwKtXr7evXdOA54DXf45eAyb7fHwOu9/1+A/CY7/fJwIu+33v4tncDoKPvc5BdVz8TeGuW/8T3ex5QlM7bGW9J3u1AQ7/t+x/ptp2Bc4EzgI/8prm+XUMtI2xfU/1HUEsbZAiw0O/x7cDtqe5XAuvzCnAhsAlo7ZvWGtjk+/1x4Aq/+Tf5nr8CeNxv+uO+aa2BjX7TA+ZL4Xq2Bf4JXAC85vvA7wNygrcr3rodQ3y/5/jmk+Btbc1XFz8TQBPfDlKCpqftdub7Ou3NfdvtNWBEOm5noAOBQcH17RpqGeH+ZcrpI+uDZ6nwTat3fIfLHmAlcJIxZjeA7/9WvtlCrW+46RUO01Ptf4DbgGrf4xbAAWNMle+xfz/tdfM9f9A3f6zvRSp1AvYCf/GdMntKRBqRxtvZGLMT+D3wObAb73ZbTXpvZ0ttbNdQywgpU4KC03nTejfsSkQaA/8AfmGMORRuVodpJo7pKSMiFwN7jDGr/Sc7zGoiPFdv1hnvN98zgEeNMR7gW7yH/KHU+3X2neMeh/eUzylAI2CUw6zptJ0jSek6ZkpQqADa+T1uC+xKUV/iIiK5eAPC340xc32TvxSR1r7nWwN7fNNDrW+46W0dpqfSWcBYEfkUeAHvKaT/AYpExCoj699Pe918zzfFW+I11vcilSqACmPMSt/jOXiDRDpv5+HAdmPMXmPMcWAucCbpvZ0ttbFdQy0jpEwJCquArr4RDXl4L1DNT3GfouYbSfAn4GNjzAN+T80HrBEI1+C91mBNv9o3imEwcNB36LgQ+IGINPN9Q/sB3vOtu4HDIjLYt6yr/dpKCWPM7caYtsaYDni317+MMT8ElgCX+mYLXmfrvbjUN7/xTZ/sG7XSEeiK96JcnftMGGO+AHaIyGm+ScOADaTxdsZ72miwiBT4+mStc9puZz+1sV1DLSO0VF5kquWLPKPxjtrZBsxIdX9i7PvZeA8H1wJrfP9G4z2X+k9gi+//5r75BXjEt67rgBK/tqYAW33/rvWbXgJ85HvNwwRd7Ezx+g/l+9FHnfD+sW8FXgYa+Kbn+x5v9T3fye/1M3zrtQm/0TZ18TMB9APKfNt6Ht5RJmm9nYF7gI2+fj2LdwRRWm1n4Hm810yO4/1m/+Pa2K6hlhHun2Y0K6WUsmXK6SOllFJR0KCglFLKpkFBKaWUTYOCUkopmwYFpZRSNg0KSimlbBoUlFJK2TQoKJUEIjJPRFb76gJcl+r+KBUvTV5TKglEpLkx5isRaYj31grnGWP2p7pfSsUqJ/IsSqko/KeIXOL7vR3ee+9oUFD1jgYFpRIkIkPx3u1ziDGmUkSW4r1Hj1L1jl5TUCpxTYGvfQGhO94SikrVSxoUlErcm0COiKwF7gNWpLg/SsVNLzQrpZSy6ZGCUkopmwYFpZRSNg0KSimlbBoUlFJK2TQoKKWUsmlQUEopZdOgoJRSyqZBQSmllO3/A7XnAX+MjiHjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from CheckPrime import is_prime\n",
    "from ModulesFactorization import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "print(\"n^2-aが素数になる確率を計算します\")\n",
    "\n",
    "a_min_str=input(\"自然数aの最小値を入力してください:\")\n",
    "a_min=PositiveInt(a_min_str)\n",
    "\n",
    "a_max_str=input(\"自然数aの最大値を入力してください:\")\n",
    "a_max=PositiveInt(a_max_str)\n",
    "\n",
    "num_str=input(\"各aに対して調べる自然数の個数Numを入力してください:\")\n",
    "num=PositiveInt(num_str)\n",
    "    \n",
    "a_list=[]\n",
    "Pa_list=[]\n",
    "max_Pa=0\n",
    "max_a=0\n",
    "\n",
    "for a in range(a_min,a_max+1):\n",
    "    if not is_square(a):\n",
    "        ini_n=round(a**0.5)+1\n",
    "        a_list.append(a)\n",
    "        count_prime=0\n",
    "        for n in range(ini_n,ini_n+num):\n",
    "            p=n**2-a\n",
    "            if is_prime(p):\n",
    "                count_prime+=1\n",
    "                #print(p)\n",
    "        Pa=count_prime/num\n",
    "        if Pa > max_Pa:\n",
    "            max_Pa=Pa\n",
    "            max_a=a\n",
    "        Pa_list.append(Pa)\n",
    "        #print(count_prime/num)\n",
    "\n",
    "        \n",
    "#print(a_list,Pa_list)\n",
    "plt.plot(a_list,Pa_list,color=\"#77ff77\",marker=\"o\",linestyle=\":\") #エメラルドグリーンでプロット\n",
    "plt.plot(max_a,max_Pa,color=\"red\",marker=\"v\",label=\"a={0},Pa={1}\".format(max_a,max_Pa))\n",
    "\n",
    "plt.xlabel(\"a\")\n",
    "plt.ylabel(\"P(a)\")\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig(\"n2a_a_min={0}_a_max={1}_Num={2}\".format(a_min,a_max,num))\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

	# coding: utf-8

	# In[29]:


	def is_prime(num_in):
	if num_in <= 1:
	return False
	else:
	bound=round(num_in**0.5)
	for num in range(2,bound+1):
	if num_in % num ==0:
	return False
	return True

	def is_prime2(num_in):
	if num_in <= 1:
	return False
	if num_in%2 ==0:
	return True
	else:
	bound=round(num_in**0.5)
	bound=(bound-1)//2
	for num in range(1,bound+1):
	if num_in % (2*num+1) ==0:
	return False
	return True

	import random

	def is_prime3(n):
	"""
	return True if num_in is probably prime by Miller-Rabin primary test.
	: param n : int (natural number)
	"""
	NumTrial = 100 # Miller-Rabin法での試行回数、素数判定を誤る確率は4^(-NumTrial)

	if n <= 0: return False
	if n == 2: return True
	if n == 1 or n & 1 == 0: return False

	d = (n - 1) >> 1
	while d & 1 == 0:
	d >>= 1

	for k in range(NumTrial):
	a = random.randint(1, n - 1)
	t = d
	y = pow(a, t, n)

	while t != n - 1 and y != 1 and y != n - 1:
	y = (y * y) % n
	t <<= 1

	if y != n - 1 and t & 1 == 0:
	return False

	return True

	# coding: utf-8

	# In[9]:


	import sys

	def result_factorization(N,factors):
	"""
	return True if non-trivial divisor of N is found
	:param N : int (natural number)
	:param factors : (list of factors of N, len(factors)=2)
	"""
	if len(factors) !=2:
	print("len(factors) != 2")
	sys.exit()

	if N != factors[0] * factors[1]:
	print("factors[0]×factors[1] != N @print_factors in ModulesFactorization")
	sys.exit()

	if factors[0] ==1:
	print("自明な約数しか見つかりませんでした。")
	return False
	else:
	print("{0}={1}×{2}".format(N,factors[0],factors[1]))
	return True


	# In[1]:


	import sys
	import math

	def is_square(N):
	"""
	return True if N is square number
	:param N : int (natural number)
	"""
	if N < 0:
	print("N is negative number @is_square in ModulesFactorization.")
	sys.exit()

	sqrt_N=round(math.sqrt(N))
	if N == sqrt_N*sqrt_N:
	return True
	else:
	return False


	# In[2]:


	import sys
	import math
	import numpy as np

	def GenerateCoprimes(N):
	"""
	retruns list of numbers which is relatively prime to N
	:param N :int (Natural number)
	"""
	if N <= 0:
	print("N <= 0 @GenerateCoprimes in ModulesFactorization")
	sys.exit()

	Coprimes = list() # Nまでの数と互いに素な数のリスト
	for num in range(1,N+1):
	if math.gcd(N,num) == 1:
	Coprimes.append(num)
	return Coprimes


	# In[3]:


	def FactorInList(N, List):
	"""
	return factor of N in PrimeList if exists
	:param N: int (Natural number)
	:param List: list of int (Candidate of factors)
	"""
	a = 1
	b = N
	for num in List:
	if N % num == 0:
	a = num
	b = N//num
	break

	return [a, b]


	# In[4]:


	import sys

	def PositiveInt(N_str):
	"""
	return positive integer N for string N_str
	:param N_str: string
	"""
	try:
	N = int(N_str)
	except ValueError:
	print("整数を入力してください。")
	sys.exit()

	if N <= 0:
	print("0以下の整数です。1以上の自然数を入力してください。")
	sys.exit()

	return N