{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e279c7a3-8b47-4d07-9ffa-80462208061d",
   "metadata": {},
   "source": [
    "# congruent conditional prime distribution\n",
    "### Hiroyuki OGAWA (Osaka Univ.)\n",
    "#### The original version of this program was written in Mathematica and made for a talk at the RIMS conference.The original version of this program was written in Mathematica and made for a talk at the RIMS conference on 2013/11/5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7971794-6f61-44ba-b4c4-f33977a835fe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import primes.primes as primes\n",
    "import matplotlib.pyplot as plt\n",
    "from math import log10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6396669-0f7d-4810-8fca-4ecb3960e416",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def chop(x, e=3):\n",
    "    return(int(100*x*10**e+1/2)/10**e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e8eac75-bd4d-4085-8df5-36715428f4ac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ccpd(X, d, m, residue=()):\n",
    "    if type(m)==int:\n",
    "        m = (m, m)\n",
    "    if len(m)<2:\n",
    "        m = m*2\n",
    "    if len(residue) < 2:\n",
    "        residue = (primes.irreducible_residues(m[0]), \n",
    "                   primes.irreducible_residues(m[1]))\n",
    "    prsm = [(p%m[0], p%m[1]) for p in primes.table((X, X+d))]\n",
    "    ccM = [[0]*m[1] for _ in range(m[0])]\n",
    "    s = prsm[0]\n",
    "    for t in prsm[1:]:\n",
    "        ccM[s[0]][t[1]] += 1\n",
    "        s = t\n",
    "    n = {k:sum(ccM[k]) for k in residue[0]}\n",
    "    return({'modulo':m, 'residue':residue, 'X':(X,d), \n",
    "            'data':{k:{j:(ccM[k][j], chop(ccM[k][j]/n[k])) \n",
    "                       for j in residue[1]} \n",
    "                    for k in residue[0]}})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb88e10b-4484-434e-bbba-9babe96a8a86",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def listTranspose(a):\n",
    "    return([list(x) for x in zip(*a)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de950b3c-d663-4e9c-aef8-0be41635859c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ccpdAll(m, Xs, residue=()):\n",
    "    if type(m)==int:\n",
    "        m = (m, m)\n",
    "    if len(m)<2:\n",
    "        m = m*2\n",
    "    if len(residue) < 2:\n",
    "        residue = (primes.irreducible_residues(m[0]), \n",
    "                   primes.irreducible_residues(m[1]))\n",
    "    wk = {X:ccpd(*X[-2:], m, residue=residue)['data'] for X in Xs}\n",
    "    return({'modulo':m, 'residue':residue, 'Xs':Xs, \n",
    "            'data':[[((i,j),listTranspose([wk[X][i][j] for X in Xs])) \n",
    "                     for j in residue[1]] \n",
    "                    for i in residue[0]]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd0abe13-6798-458a-95d8-c79cd1c2731b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ccpdPlot(m, Xs, residue=(), size=10):\n",
    "    if type(m)==int:\n",
    "        m = (m, m)\n",
    "    if len(m)<2:\n",
    "        m = m*2\n",
    "    if len(residue) < 2:\n",
    "        residue = (primes.irreducible_residues(m[0]), \n",
    "                   primes.irreducible_residues(m[1]))\n",
    "    wk = ccpdAll(m, Xs, residue=residue)['data']\n",
    "    idx = [log10(max(1,X[-2])+X[-1]/2) for X in Xs]\n",
    "    #\n",
    "    for ds in wk:\n",
    "        c = 100/len(ds)\n",
    "        plt.plot(idx, [c]*len(idx), c='red')\n",
    "        for ij,data in ds:\n",
    "            plt.scatter(idx, data[1], s=size, \n",
    "                        label='{} ({})'.format(ij[1],m[1]))\n",
    "        plt.yticks([c-5,c-1,c,c+1,c+5], \n",
    "                   ['-5%','-1%','mean','+1%','+5%'])\n",
    "        plt.legend(title='{} ({}) →'.format(ij[0],m[0]),\n",
    "                   loc=\"upper left\", bbox_to_anchor=(1.02, 1.0,), \n",
    "                   borderaxespad=0, handletextpad = 0)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7431c7b3-8d90-4987-8cf9-7e1260f501c5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Xs = sum([[(4**i+j*10**5, 10**5) for j in range(10)] for i in [0,8,9]]\n",
    "         +[[(4**i+j*10**6, 10**6) for j in range(10)] for i in [0,10,11]]\n",
    "         +[[(4**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(12, 14))]\n",
    "         +[[(2**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(28, 51))], [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b75aa72-6f44-48c9-9b4e-c285fef1e83c",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xs8 = sum([[(i, 2**i+2*j*10**8, 10**8) for j in range(5)] for i in range(28, 51)], [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a3a0453-5cdd-48f8-87d0-46b8825a8a8f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Xs = sum([[(4**i+j*10**5, 10**5) for j in range(10)] for i in [0,8,9]]\n",
    "#         +[[(4**i+j*10**6, 10**6) for j in range(10)] for i in [0,10,11]]\n",
    "#         #+[[(4**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(12, 14))]\n",
    "#         #+[[(2**i+j*10**7, 10**7) for j in range(10)] for i in [0]+list(range(28, 51))]\n",
    "#         , [])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27444e15-c4d9-4bdf-92bf-f7efd522acf9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(4, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0894d43e-0d13-4252-a9e9-928e5c1db232",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(6, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61eef572-9f46-4f4d-a818-058354352b99",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(5, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b40777-0b2c-48cd-9ca7-ab8b9d562e51",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(7, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f247d4d3-2842-47fa-bc3e-161cb0ad8485",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(8, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18e22134-5893-406c-9755-5bd8937688e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(9, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e8a3cb-71ed-48c0-8829-2eef3c323354",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot(11, Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f30a0990-a309-40cd-b598-8e90192f58e9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((4,6), Xs)\n",
    "ccpdPlot((6,4), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a865228-6e8e-45b2-adf7-6577a05b092e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((4,5), Xs)\n",
    "ccpdPlot((5,4), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1fd0d98-5c02-40f9-bf79-ed43c2699e47",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((6,5), Xs)\n",
    "ccpdPlot((5,6), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4fa7809-f3b7-4b44-aac1-69791fc09da6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((5,7), Xs)\n",
    "ccpdPlot((7,5), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74206483-780d-4d4c-9898-3c8ffc432b45",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((5,7), Xs8)\n",
    "ccpdPlot((7,5), Xs8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d324af5-bd79-4fce-b595-827f946abe1a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((4,8), Xs)\n",
    "ccpdPlot((8,4), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb85adc3-b49f-4e51-b8a7-4591239ee8aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((6,8), Xs)\n",
    "ccpdPlot((8,6), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d80c2ecf-0a34-4d39-bdbb-400b12aee52c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ccpdPlot((6,7), Xs)\n",
    "ccpdPlot((7,6), Xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cd7cc06-f7a0-41ac-b261-466566b077b5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01b73057-3dc4-4365-8bc5-9a7862ba5877",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (HomeBrew)",
   "language": "python",
   "name": "python3-hb"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}