four_ducks_Monte.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "author: 赵杭天 @ SINAP\n",
    "date: 2019.8.86\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 根据你的计算机性能时间修改采样数,这决定了计算时间(sample_num)\n",
    "sample_num = 10000*16*4\n",
    "xx = np.array([random.random()-0.5 for i in range(sample_num)])\n",
    "yy = np.array([random.choice([1, -1])*np.sqrt(0.5**2-xx[i]**2)*random.random() for i in range(len(xx))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.2008448897813685, 0.03962637862371878)\n",
      "(-0.34010471319945634, -0.20119591994657227)\n",
      "(-0.29679665764259955, -0.0968709610924799)\n",
      "(0.26836952735581354, 0.3606777356207143)\n",
      "(0.33039570717006417, -0.018156378167732166)\n",
      "(0.16342082402698832, -0.40061986749355666)\n",
      "(0.47242536365121623, -0.014152163902939717)\n",
      "(0.3808843690325415, 0.1360950737641591)\n",
      "(0.26678817796624776, 0.09584330358204363)\n",
      "(-0.4680732558488545, -0.04141602740924838)\n",
      "(0.22507176615613678, 0.05328757454270317)\n",
      "(-0.130169666192666, -0.40441729170684665)\n",
      "(0.06333268541990644, 0.029594204785379472)\n",
      "(-0.16113789559275216, -0.24295417559694196)\n",
      "(0.3838365492117176, 0.31486337882682985)\n",
      "(-0.39362056105191523, 0.07259820049110674)\n"
     ]
    }
   ],
   "source": [
    "# 打印前16组点\n",
    "for i in range(0,16):\n",
    "    print((xx[i], yy[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_exist_solve(xx, yy):\n",
    "    theta = np.arange(0.0,np.pi,0.01)\n",
    "    p1 = np.array([xx[0],yy[0]])\n",
    "    p2 = np.array([xx[1],yy[1]])\n",
    "    p3 = np.array([xx[2],yy[2]])\n",
    "    p4 = np.array([xx[3],yy[3]])\n",
    "    for t in theta:\n",
    "        p0 = np.array([0.5*np.cos(t), 0.5*np.sin(t)])\n",
    "        # 向量法判定点相对直线的位置\n",
    "        if (np.cross(p1,p0)>=0. and np.cross(p2,p0)>=0. and  np.cross(p3,p0)>=0. and np.cross(p4,p0)>=0.) or (np.cross(p1,p0)<=0. and np.cross(p2,p0)<=0. and  np.cross(p3,p0)<=0. and np.cross(p4,p0)<=0.):\n",
    "            return True\n",
    "    return False\n",
    "\n",
    "def get_exist_solve(xx, yy):\n",
    "    theta = np.arange(0.0,np.pi,0.01)\n",
    "    p1 = np.array([xx[0],yy[0]])\n",
    "    p2 = np.array([xx[1],yy[1]])\n",
    "    p3 = np.array([xx[2],yy[2]])\n",
    "    p4 = np.array([xx[3],yy[3]])\n",
    "    for t in theta:\n",
    "        p0 = np.array([0.5*np.cos(t), 0.5*np.sin(t)])\n",
    "        # 向量法判定点相对直线的位置\n",
    "        if (np.cross(p1,p0)>=0. and np.cross(p2,p0)>=0. and  np.cross(p3,p0)>=0. and np.cross(p4,p0)>=0.) or (np.cross(p1,p0)<=0. and np.cross(p2,p0)<=0. and  np.cross(p3,p0)<=0. and np.cross(p4,p0)<=0.):\n",
    "                print('found solution')\n",
    "                return p0\n",
    "    print(\"can't find solution\")\n",
    "    return 0,0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "can't find solution\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEWCAYAAAC9hIj8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW9//HXZ7LvITshIQskQUDWsKm4IRasorbuetWWSq2/Vq221lZ729vWtmpr7+212lLrdaPuS1GpS1VcICxhlS0sgZCQkH1fJst8f39ksAETEpKZOZPM5/l45JHMzMk5H0Lmne/5nvP9fsUYg1JKuYPN6gKUUiOXBoxSym00YJRSbqMBo5RyGw0YpZTbaMAopdxGA0aNOCLycxF5zuo6lAaMTxORs0RkrYjUi0iNiKwRkVlD3OfNIvLZCc89JSK/Glq1XzrOUyLSLiJNztrfF5EJg9jPIRG5wJW1qX/TgPFRIhIJvAX8LxADjAH+C7BbWVdvRMS/j5ceMsaEAylABfCUx4pSA6IB47uyAYwxzxtjuowxrcaY94wx249tICK3iMhuEWkUkV0iMsP5/L0icqDH85c7nz8N+DMwz9myqBORZcD1wD3O5950bpssIq+KSKWIHBSR23sc9+ci8oqIPCciDcDNJ/uHGGNagL8Dk3t7XUSWiMhOZz2rnXUiIs8CY4E3nbXdM7gfpeqLBozv2gt0icjTIrJYREb1fFFErgR+DtwIRAJLgGrnyweA+UAU3a2e50RktDFmN3ArkGeMCTfGRBtjlgMrcLY2jDGXiIgNeBPYRnfLaQFwp4h8pUcJlwKvANHO7++TiITTHWJbenktG3geuBOIB1bRHSiBxpj/AA4Dlzhre6j/H5s6FRowPsoY0wCcBRjgr0CliKwUkUTnJt+iOxQ2mm77jTFFzu992RhTaoxxGGNeBPYBs0/h8LOAeGPML4wx7caYQmcN1/TYJs8Y84bzGK197OcHIlIH7AfC6b2lczXwtjHmfWNMB/A7IAQ44xTqVYPU17mt8gHOFsfNAM4O0ueA/wauBVLpbql8iYjcCNwFpDufCgfiTuHQaUCyMxyO8QM+7fG4eAD7+Z0x5v5+tkkGio49MMY4RKSY7paTcjMNGAWAMWaPiDwFfNv5VDEw7sTtRCSN7tbGArpbGV0ishWQY7vqbfcnPC4GDhpjsk5W0imUfzKlwOnHHoiI0B2eR1x8HNULPUXyUSIyQUTuFpEU5+NUulsu65ybPEH3KchM6TbeGS5hdL8pK53f9w2O71wtB1JEJPCE5zJ7PN4ANIjIj0QkRET8RGTyUC+R9+El4KsiskBEAoC76b5StraP2pQLacD4rkZgDrBeRJrpDpYddL8BMca8DDxA99WZRuANIMYYswv4PZBH95vzdGBNj/1+COwEjopIlfO5vwETnVdx3jDGdAGXANOAg0AV3YEW5ep/pDGmALiB7svxVc7jXmKMaXdu8hvgfmdtP3D18X2d6IRTSil30RaMUsptNGCUUm6jAaOUchsNGKWU24y4+2Di4uJMenq61WUoNaJt2rSpyhgT3992Iy5g0tPTyc/Pt7oMpUY0ESnqfys9RVJKuZEGjFLKbTRglFJuM+L6YJQaLjo6OigpKaGtrc3qUvoUHBxMSkoKAQEBg/p+DRilLFJSUkJERATp6el0D/L2LsYYqqurKSkpISMjY1D70FMkpSzS1tZGbGysV4YLgIgQGxs7pBaWBoxSFvLWcDlmqPXpKZJym44uBxWNdqoa7TTZO2ls66TJ3klTWwetHQ4czpH8xhiMAZtNCA30IyzQn5BAP8KC/IgMDiA+Ioj4iCBCA/XXdbjR/zE1aMYYaprbOVDZTGFlEwcqmzhY1cLRhlaO1tupbrbjytlAwoP8SYgIYnR0MOmxYWTE/ftjbEwo/n7aIPc2GjBqwI7Wt7GtpI7tJXVsL6lnx5F6als6vng9yN9GemwYo6ODmZwcRWJkMElRwcSHBxEZEkBYkB8RQQGEB/sTEuCHCN0fCCLQ5TC0tnfR0tFFi72T5vYu6ls7qGy0U9HY5vxs50htK29tL6O+9fhjTxgdyeljIjl9TBSTx0SRkxihoWMxDRjVp9K6VtbsryLvQDV5hdWU1Xd39vnZhOzECC6cmER2UgTj4sMYFx/OmOgQbLbBn7MH+EFwgB+j+t8UgNrmdgqrultPBUcb+fxIPW9sKeW5dYcBCAv0Y1ZGDHMzY5mXGcuk5MhhHTjl5SsoLLwPu/0wQUFjycx8gMTE6we9v40bN7J06VI2bNhAV1cXs2fP5sUXX2Ty5F6XlxoUDRj1hY4uB+sLa3hv11E+2VvJoeoWAGLDApk7LpbctFFMSYlmUnIkwQF+FlcLo8ICmRkWyMy0f0eSw2Eoqmlhe0kdGw/VsK6wht/+cw/QfYo1PyuOhRMTOS8ngVFhgX3t2uuUl6+goGAZDkf3/4ndXkRBwTKAQYfMrFmzWLJkCffffz+tra3ccMMNLg0XGIFTZubm5hod7Dhwre1dfFRQwXs7j/Lhngoa2joJCfDjjHGxnDE+jjPHx5KdEDGklonVKhrb2HCwhjX7q/hgdwUVjXZsArnpMVw4MZGLpySTFBXs8bp2797NaaedNqBt8/LSsdu/PL4wKCiNefMODbqG9vZ2Zs2aRXBwMGvXrsXP78t/OHqrU0Q2GWNy+9u/tmB8kMNhWH+whtc2l7Dq8zKa27sYFRrAhZOS+MqkJM4aH0dIoPUtFFdJiAjm4inJXDwlGYfDsKO0nn/tKuf93RX86u3dPLBqN2eOi+Oy6WNYNDmJ8CDve1vY7YdP6fmBqqmpoampiY6ODtra2ggLCxvS/k6kLRgfUlbfyvMbinltcwklta2EBfpx0emjuXz6GGZnxAzr/onBOljVzBtbjvDG1iMUVbcQHGBj8eTR3DA3jRljo916n4o3tGCWLFnCNddcw8GDBykrK+PRRx8dUJ3aglFA96Xk/KJanlpziHd2HsVhDGeNj+MHF+bwlUlJI6qlMhgZcWF8f2E2d16QxebDdby+pYQ3tpTy+pYjTEqO5MZ5aSyZOsbyn1Nm5gPH9cEA2GyhZGY+MOh9PvPMM/j7+3PdddfR1dXFGWecwYcffsj555/vipIBbcGMWJ1dDlZuK+Vvnx1kZ2kDkcH+XDt7LDfMTSM1JtTq8rxas72T17cc4dm8IgrKG4kKCeDGeWl848wMYlzYMXwqLRhw/VWkgRpKC0YDZoRp73Tw2uYSHlt9gMM1LWQnhnPzGRlcNj1Z74Q9RcYYNhys4ck1B3l3ZzkhAX5cN2cst8zPdEmn8KkGjFX0FEnR3ungxfxi/rz6AEfqWpmSEsVPL87lgtMSvH68i7cSEeZkxjInM5Z95Y08vvoAT609xLN5RVw1K4Xbz88iIdLzV5+GEw2YYc4Ywzs7jvLgO3s4VN3CzLRRPHD5ZM7JjtdgcaGsxAgeuXoa31+YzeMfH+CFDcW8uukIt8zP4JazM4kIHtx8KcYYr/5/GuoZjp4iDWObimr59ardbCqqJTsxnB8vPo1zczRYPOFQVTO/e6+At7aXERsWyPfOH88Nc9NO6UrcwYMHiYiI8NopG47NB9PY2Pil+WC0D2YEq2y088Dbu3hjaykJEUHctTCbK2am+ORlZqttL6njt//cw9oD1UxIiuBXl00mNz1mQN87nGe004AZgRwOw983HObBd/Zg73Bw6zmZ3HruOO28tZgxhnd3lvOLN3dSWt/GlTNTuHfxBGLDg6wuzW20k3eE2V3WwI9f+5ytxXWcMS6WX142mXHx4VaXpejuDF40OYmzs+P44wf7eeLTQt7bVc5/XjyRr80Y45WnP56iLRgv1+Uw/OWTA/zh/b1EBgdw/8Wncdk03/6l9Xb7Kxr58Wufs/FQLQsnJvLry08nPmJktWYG2oKx9KRdRBaJSIGI7BeRe0+y3RUiYkSk33/QSFJU3czVf8njoXcKWDgxkffvOofLp6douHi58QkRvLBsHvd/9TQ+3lvJhX/4mLe3l1ldliUsCxgR8QP+BCwGJgLXisjEXraLAG4H1nu2QusYY3hpYzGL/+dTCsob+cPVU/nTdTNcehepci8/m/Ct+Zmsuv0sUmNC+X9/38zdL22jpb3T6tI8ysoWzGxgvzGm0BjTDrwAXNrLdr8EHgK8t6vdhVrbu/jhK9u559XtTEuN5t07z9ZWyzA2PiGC175zBrefP57XtpRw6aNr2FfeaHVZHmNlwIwBins8LnE+9wURmQ6kGmPeOtmORGSZiOSLSH5lZaXrK/WQwsomLn9sDa9uLuGOBVk8u3QOydEhVpelhsjfz8ZdF+bwzDdnU9PczpJH1/Da5hKry/IIKwOmtz/JX/Q4i4gN+ANwd387MsYsN8bkGmNy4+PjXVii57y78yhLHl1DeUMbT31jNt9fmI3fMJ7kSX3Z/Kx4Vt0xn9NTorjrpW387B876OxyWF2WW1kZMCVAao/HKUBpj8cRwGRgtYgcAuYCK0daR68xhsdW7+fbz25iXEI4b90+n3Oyh2dIqv4lRgbz92/N4Zb5GTydV8Q3ntp43OTlI42VAbMRyBKRDBEJBK4BVh570RhTb4yJM8akG2PSgXXAEmPMiLkG3d7p4J5XtvPQOwVcMjWZF5fNZYyeEo14/n427vvqRB76+hTWFVZz+WNrOFTVbHVZbmFZwBhjOoHvAu8Cu4GXjDE7ReQXIrLEqro8pb6lgxufXM/Lm0q4fUEWf7xmmldMpK0856pZqTy3dA61ze1c+qc1bCqqtbokl9Mb7SxQ0dDGf/xtAwermnnwitO5fHqK1SUpCx2ubuHGJ9dT3mDn8RtmcG5OgtUl9WtY3Gjni4prWrjyL3kU17bwf9+YpeGiGBsbysu3nkFGXBjfejqfldtK+/+mYUIDxoP2VzRx1V/yqG1u59mlczhzfJzVJSkvER8RxAvfnsuMsaO444UtPLfuyxN8D0caMB6yr7yRq/+SR0eX4cVvzztusTClACKDA3hm6WzOy0ng/jd28Pf1Q1uSxBtowHjAwapmrntiPTab8NK353La6EirS1JeKjjAj8dvmMF5OfH85PXPeWljcf/f5MU0YNysuKaF6/+6ji6H4e/fmkOmTrGg+hHk78fjN8xkflYcP3pt+7C+61cDxo2O1rdx/RPrabJ38uzS2WQlRlhdkupDefkK8vLSWb3aRl5eOuXlKyytJzjAj7/emMu8zFh+8PI23t151NJ6BksDxk0a2zq4+f82UN1k5+lvzmZScpTVJak+HFtYvnvlRPPFwvLeEDJP3JTLlJRobn9+C5uKaiytZzA0YNygo8vBbSs2s6+iicdumMn0sdqh680KC+87bsVEAIejhcLC+yyq6N9CA/352025jI4KZunT+RyobLK6pFOiAeNixhjuf30Hn+6r4tfO5UOUd3PXwvKuEhsexNPfnI2/TbjpyQ1UNA6fmUs0YFzssdUHeDG/mO+dP56rZ421uhw1AEFBvf8/9fW8FdJiw3jy5llUN7Xz7Wc3Ye/ssrqkAdGAcaGPCir43XsFLJmazF0Ls60uRw1QZuYD2GzHr9c91IXl3WFKSjS/v2oqWw7X8fOVu6wuZ0A0YFzkcHULd76wlQlJkTz49Sk6A90wkph4PTk5ywkKSgOEoKA0cnKWe2Rh+VN10emjue3ccTy/4fCwuBFPly1xgdb2Lm59bhPGGP58wwxCAnVU9HCTmHi9VwZKb+6+MIedpQ38bOUOcpLCmZk2sIXerKAtGBf46T92sPtoA/9zzXTSYsOsLkeNcH424Y/XTGd0VAjf+/sW6lu8d8IqDZghenNbKa9sKuG7543nvAneP8xejQxRoQH877XTqWi08+PXtw95kXp30YAZgiN1rdz3+udMHxvNHQuyrC5H+ZipqdHcfWEOqz4/yoteOmZJA2aQuhyGu17cSpfD8N9XT9OF55Ulvn12JmeOj+W/3tzF/grvuwlP3xWD9MSnhaw/WMPPl0zSfhdlGZtNeOSqaQQH2PjhK9vocnjXqZIGzCAcqmrmkff3cuHERK6YqTPSKWslRgbzs0smseVwHU+vPWR1OcfRgDlFxhh+8vrnBPrb+OVlk/V+F+UVLp2WzLk58Tz8bgHFNS39f4OHaMCcopc3lbD2QDX3Lp5AYmSw1eUoBYCI8MDlp2MT+Mnrn3vNVSUNmFNQ1WTngbd3Myt9FNfqOCPlZcZEh/CjxRP4dF8Vb24vs7ocQAPmlPz+vb002zv5zddOx6bLuiovdP2cNCYlR/KbVbtpbbd+QKQGzADtLmvgxY2H+Y95aYxP0JnplHfyswk/u2QSZfVt/PnjA1aXowEzEMYYfvX2LiJDAvSGOuX1ZmfEcPGU0fz54wOU1Frb4asBMwD/2l3Bmv3V3Lkgi+jQQKvLUapfP77oNETgoXcKLK1DA6YfXQ7Dg+/sITM+jOvnplldjlIDMiY6hG+emcGb20vZc7TBsjo0YPrx1vZS9lc0cdfCbAJ0OIAaRpadnUl4oD9/eH+vZTXoO+YkuhyGP36wj5zECC6aPNrqcpQ6JdGhgSydn8G7O8v5vKTekho0YE7izW2lHKhs5o4LsvSytBqWvnlWBtGhATzyvjV9MZYGjIgsEpECEdkvIvf28vpdIrJLRLaLyAci4rFOEIfD8McP9zEhKYJFk5I8dVilXCoyOIBb5mfyUUElu8s83xdjWcCIiB/wJ2AxMBG4VkQmnrDZFiDXGDMFeAV4yFP1fbingsLKZm47b7y2XtSwdv2csYQG+vHXTws9fmwrWzCzgf3GmEJjTDvwAnBpzw2MMR8ZY45dyF8HeGzo8t8+O0hyVDAXTdbWixreokMDuSo3lTe3lXK03rNrKlkZMGOAntNwlTif68tS4J9urchpZ2k9eYXV3HRGuk4kpUaEpWdl0OUwPOXh6RysfPf0dt7R6xBQEbkByAUe7uP1ZSKSLyL5lZWVQy7sb58dJDTQj2tm64BGNTKkxoSyePJo/r6+yKNjlKwMmBIgtcfjFKD0xI1E5ALgPmCJMcbe246MMcuNMbnGmNz4+KEt1Vrb3M5b28q4YmYKUSEBQ9qXUt7khrlpNLR1supzz420tjJgNgJZIpIhIoHANcDKnhuIyHTgL3SHS4Uninpj6xHauxxcq60XNcLMzYwhIy6MFzZ6bsE2ywLGGNMJfBd4F9gNvGSM2SkivxCRJc7NHgbCgZdFZKuIrOxjd66qiRc3FjMlJYrTRke681BKeZyIcPWsVDYeqvXYBOGW9mAaY1YZY7KNMeOMMQ84n/tPY8xK59cXGGMSjTHTnB9LTr7HodlxpIE9Rxu5Mje1/42VGoa+PiMFf5vwoodaMXqJpIeX8osJ8rexZGqy1aUo5RbxEUEsOC2BN7aWemQFAg0Ypy6HYdXnZVwwMVE7d9WIdvGUZCob7Ww8VOP2Y2nAOG04WEN1cztfPV0HNaqRbcFpCQQH2HjbA/P2asA4vbOjjOAAG+fmDO0yt1LeLjTQnwUTEvnnjjI6uxxuPZYGDN0DG/+54yjnZicQGuhvdTlKud1Xp4ymqqmdDQfde5qkAQNsKa6jotHO4tN13JHyDefmxBPoZ+OjAvfeXqYBA3yytxKbwLnZCVaXopRHhAb6Mzsjho/3Dn1ozclowACf7qtkSko0UaF69Uj5jnOy49lb3kRpXavbjuHzAVPf2sHW4jrOzoqzuhSlPOrYBQ13tmJ8PmDyDlThMDA/W68eKd8yPiGc5KhgPt2nAeM2aw9UExrox7TUaKtLUcqjRITZGTFsPFSLMe65q9fnA2bz4VqmpUbrkiTKJ81Mj6Gy0U5xjXv6YXz6XdXa3sXuskZmjB1ldSlKWSI3rft3f9Nh99wP49MBs72kji6HYfpYPT1Svik7MYKIIH/yD9W6Zf8+HTBbiusAtP9F+Sw/mzBtbDRbDte5Zf8+HTCfH6knNSaE2PAgq0tRyjKnjY5kf2WTW8Yl+XTA7D3aSE5ihNVlKGWpCUkRtHc6OFTd7PJ9+2zAtHc6OFjVTJYGjPJxOUnd74E9Rxtdvm+fDZhD1c10OgzZieFWl6KUpcYnhONnEwo0YFxnb3n3DzMrQVswyrcF+fuRFhvKvnLXTwTuswFzqKr7fDMzPsziSpSy3tiYUI64YdCjzwbMkbo2YsICdYIppYCUUSEU17b0v+Ep8tmAKa1rJTk62OoylPIKKaNCqWvpoLGtw6X79dmAKatvJTkqxOoylPIKKaO63wuuPk3yyYAxxlBa10ZytAaMUgCjnX9sy+rbXLpfnwyYtg4HTfZO4iP0Dl6lAEY5Z3Osa2l36X59MmDqW7vPM6N1ikylAIgJCwSgpln7YIbsWMDoCo5KdYsMDsAm2oJxCQ0YpY5nswlRIQHUasAMnQaMUl8WERxAY1unS/fpkwHT2tEFQGign8WVKOU9Av1tdHa5dm5eSwNGRBaJSIGI7BeRe3t5PUhEXnS+vl5E0l1x3C5H97wX/jafzFelehXgZ8Pe6do5Yfp8h4nIKle9ofvYvx/wJ2AxMBG4VkQmnrDZUqDWGDMe+APwoCuOfSyl/Wziit0pNSIE+tvocPGkUyf7E/4U8J6I3Cci7uismA3sN8YUGmPagReAS0/Y5lLgaefXrwALRGTIqdDp6A4Yfz8NGKWO2VZc5/JF2Poc6WeMeUlE3gb+E8gXkWcBR4/XHxnisccAxT0elwBz+trGGNMpIvVALFDVcyMRWQYsAxg7dmy/Bz4WMNqCUcq9+htK3AE0A0FABD0CxgV6e3ef2MM0kG0wxiwHlgPk5ub220vl52wEOVw/BalSw1Zu2iiCAlzbL9lnwIjIIuARYCUwwxjj6rHcJUBqj8cpQGkf25SIiD8QBQx5AZdg5w/R3tk11F0pNWJ0OAyhLr7wcbIWzH3AlcaYnS494r9tBLJEJAM4AlwDXHfCNiuBm4A84ArgQ+OCNS6D/LsvT7d1aBNGqWPaOx0Eurhf8mR9MPNdeqQv779TRL4LvAv4AU8aY3aKyC+AfGPMSuBvwLMisp/ulss1rjj2sRZMW4e2YJQ6psneQUSwa6eQtXQ6N2PMKmDVCc/9Z4+v24ArXX3c4IBjLRgNGKWOaWzrJCLYtZHgk3eaRQZ3X3U/NmRAKV9njKGxrZPwIA2YIYsJPzY03bUDu9ylvHwFeXnprF5tIy8vnfLyFVaXpEaY1o4uuhyGiGDX3vLmkzNexzrnvqgeBgFTXr6CgoJlOBzdF/Hs9iIKCpYBkJh4vZWlqRGktsU9A4B9sgUTHOBHWKAf1U3eHzCFhfd9ES7HOBwtFBbeZ1FFaiQ66pwqc3SUayfC98mAAYgND6KqyW51Gf2y2w+f0vNKDcaxgEnSgHGN0VHBlLphoSlXCwrqfehDX88rNRhl9d3vBW3BuEhabChFNa5faMrVMjMfwGYLPe45my2UzMwHLKpIjURH69sIDrBpH4yrjI0JpbLRTmu7d98Lk5h4PTk5ywkKSgOEoKA0cnKWawevcqni2haSo0NwwWQFx/HJq0gAY2O716Qurm0hO9G1dy+6WmLi9Rooyq32VzQxPj7c5fv16RYMwMGqZosrUcpaHV0OiqpbGJ+gAeMy4xPCEYGCo41Wl6KUpYqqm+l0GA0YVwoP8ictJpTdZQ1Wl6KUpfZXNAFowLjaxORIdg0iYPTWfTWS7CprxCYaMC53WlIkRdUtNNkHvhbMsVv37fYiwHxx676GjBquthXXkZ0YQWig66/5+HTATEyOBGDnkfoBf4/euq9GEmMM20rqmJoS7Zb9+3TAzBg7CoD8otoBf4/euq9GksM1LdS1dDA1VQPG5UaFBZKVEM7GQwOf5ldv3VcjydbiOgCmpka5Zf8+HTAAszJi2HSoli7HwKb61Vv31UiyrrCaiCB/ctx0s6kGTPooGu2d7Dk6sKtJeuu+GknW7K9m7rhY/P3cEwU+O1TgmDkZsQDkHahmUvLAmol6674aCQ5Xt3C4poVvnpnutmP4fAsmOTqE7MRwPiqosLoUpTxqzYHuBVLPyopz2zF8PmAAzstJYMPBmlO6H0ap4e6zfVUkRgYxzg2DHI/RgAHOm5BAR5fhs31V/W+s1AjQ1tHF6oIKzp+Q4PIpGnrSgAFmpo0iItifj/boaZLyDWsPVNHc3sVXJiW59TgaMECAn41zcxJ4f3c5HV26nKwa+d7ZcZSIIH/OGOe+/hfQgPnCJVNGU9Pczpr9epqkRrbOLgf/2l3BeRMSCPR3bwRowDidkxNPRLA/K7eVWl2KUm61rrCGmuZ2Fk127+kRaMB8Icjfj0WTknhvZ7muWa1GtFc2FRMZ7M/5ExLcfiwNmB6WTEumyd7Jh9rZq0aohrYO3tl5lCXTkgkO8HP78TRgejhjXBzJUcE8v0FHRquRadX2Mto6HFwxM9Ujx7MkYEQkRkTeF5F9zs+jetlmmojkichOEdkuIle7uy4/m3D1rLF8uq+KomqdDFyNPK9sKmF8QjhTU9wzevpEVrVg7gU+MMZkAR84H5+oBbjRGDMJWAT8t4i4Z9KKHq6elYqfTXh+Q7G7D6WUR+0srSe/qJarclPcenNdT1YFzKXA086vnwYuO3EDY8xeY8w+59elQAUQ7+7CkqKCWTAhgZfzi2nv1Hti1Mjx1JpDhAT4cXWu5+YusipgEo0xZQDOzyftzhaR2UAgcMADtXH93DSqm9t5Uy9ZqxGiusnOP7aV8rUZY4gKde3ysCfjtoARkX+JyI5ePi49xf2MBp4FvmGM6bVJISLLRCRfRPIrKyuHXPvZWXHkJEbwl08OYMzAJqJSyps9v+Ew7Z0OvuHGqRl647aAMcZcYIyZ3MvHP4ByZ3AcC5BerwuLSCTwNnC/MWbdSY613BiTa4zJjY8f+lmUiPDtczLZW96k0zioYa+to4tn8oqYnxXH+ATPLpNs1SnSSuAm59c3Af84cQMRCQReB54xxrzswdoAuGRqMslRwfz540JPH1opl3pxYzEVjXa+c+44jx/bqoD5LbBQRPYBC52PEZFcEXnCuc1VwNnAzSKy1fkxzVMFBvjZWDo/kw0Ha9hUNPBJwZXyJvbOLh5ffYDZ6THMy4z1+PEtCRhjTLUxZoExJsv5ucb5fL7Y7VDVAAARMklEQVQx5lvOr58zxgQYY6b1+NjqyTqvmZVKXHggD79boH0xalh6Ob+Eow1t3L4gy2OXpnvSO3lPIizIn/933njWFdawZn+11eUodUqOtV5mjI3mzPGeb72ABky/rpszljHRITz87h5txahh5Zm1RRypa+X7C7Mtab2ABky/gvz9uOOCLLaV1PPuznKry1FqQGqb2/nfD/dxTnY887Pcfn9qnzRgBuBr08cwPiGc3/xzt07loIaFP364jyZ7Jz+56DRL69CAGQB/Pxs/u2QiRdUt/O2zg1aXo9RJHaxq5tm8Iq6elUpOkmfvezmRBswAzc+KZ9GkJB79cD+lda1Wl6NUr4wx/GzlToL8bXx/YbbV5WjAnIr7vnoaDmP49ardVpeiVK/e2l7GJ3srufvCHBIigq0uRwPmVKTGhHLbueN5a3uZDiFQXqe+tYP/enMXp4+J4qYz0q0uB9CAOWW3nptJVkI4P3ntcxraOqwuR6kvPPTOHmqa7fz68tPxs1lzWfpEGjCnKMjfj4evnEp5Qxu/0VMl5SXyDlTz9w2HuemMdE730Gx1A6EBMwjTUqO5ZX4mz28o1uVmleUa2jq4+6WtpMeG8cOv5FhdznE0YAbp+wuzyYwL455XtlHX0m51OcqH/fwfOylvtPPIVVMJDfS3upzjaMAMUnCAH3+4ehoVjXZ+9Op2HUagLPH29jJe23KE7543nuljvzR3vuU0YIZgamo0P1o0gXd3lvPcuiKry1E+primhZ+8/jlTU6L47vnjrS6nVxowQ7T0rAzOzYnnl2/vZldpg9XlKB/R1tHFd1ZswmEMf7x2OgF+3vlW9s6qhhGbTfjdlVOJDgngthWbqG/RS9fK/f7rzZ3sONLAI1dNIy02zOpy+qQB4wJx4UE8dv0MjtS18r0XttDl0P4Y5T4v5xfz/IZibjt3HAsnJlpdzklpwLhIbnoMv7h0Mp/sreTBd/ZYXY4aoTYV1XLfGzs4Y1wsd3nBWKP+eNc1rWHu2tlj2V3WwPJPCpmQFMHXZqRYXZIaQQ5Xt3DLM/kkRwXz6HUz8PfSfpeevL/CYeanF09kbmYM9776OWsP6E14yjXqWzr4xlMbcBjDkzfPIiYs0OqSBkQDxsUC/Gz8+YaZpMWGsuyZTewsrbe6JDXM2Tu7uPW5TRyuaeEvN8wkMz7c6pIGTAPGDaJDA3n6m7OJCPbn5v/bSHFNi9UlqWGqs8vBHc9vJa+wmge/PoU5Fiw9MhQaMG6SHB3CM9+cTXungxuf3EBlo93qktQw43AY7nl1O+/sPMpPL544LPv0NGDcKCsxgr/dlMvR+jau/es6DRl1nPLyFeTlpbN6tY28vHTKy1d88Zoxhp+/uZPXNh/hroXZLD0rw8JKB08Dxs1y02N48uZZHKlt5ToNGeVUXr6CgoJl2O1FgMFuL6KgYBnl5SswxvDA27t5Jq+IZWdn8j0vHQYwEBowHjBvXCxP3jyL4toWrvvrOqqaNGR8XWHhfTgcx/fNORwt7D9wP/e/sYMnPjvITfPS+PHiCZataeQKGjAe0jNkrvpLHiW12vHry+z2w196rsth40/5X2PF+sPces44fr5k0rAOF9CA8agzxsXxzDfnUNlo54rH8yg42mh1ScoiQUFjj3vc4fDn8W33sLb0fO5emM2PFuUM+3ABDRiPm50Rw8u3zsNhDFf+eS35h2qsLklZIDPzAWy2UACa2iN4eOOvyC8/izvPaeV7Fi1U7w4aMBaYkBTJq985g7jwIK5/Yj2rPi+zuiTlYYmJ15OTs5y6zpn8at3DFNbn8MuLmrlz8RVWl+ZSGjAWSY0J5eVb5zEpOZLbVmzmkff34tBR2D6lpO0ifpH3G9rI4PllZ/MfZ19ldUkuZ0nAiEiMiLwvIvucn/uc609EIkXkiIg86skaPSE2PIjnl83lipkp/PGDfXxnxSaa7Z1Wl6XczBjDc+uKuHb5OiKD/Xn9tjOZlR5jdVluYVUL5l7gA2NMFvCB83Fffgl87JGqLBDk78fDV0zhpxdP5P1d5Xz98bUUVjZZXZZyk7aOLn7w8nbuf2MH88bF8vptZ5IR570TRg2VVQFzKfC08+ungct620hEZgKJwHseqssSIsLSszJ46huzOdrQxiX/+xn/2HrE6rKUix2ubuFrj63l1c0l3LEgiydvnsWoYTIqerCsCphEY0wZgPNzwokbiIgN+D3ww/52JiLLRCRfRPIrKytdXqynnJ0dzz/vmM/E5EjueGEr97yyjZZ2PWUa7owxvLa5hIv++CkltS08eXMu31+Y7TWrL7qT2yacEpF/AUm9vHTfAHdxG7DKGFPc3yU7Y8xyYDlAbm7usO4pHR0VwvO3zOV/PtjHox/tZ/PhOh65aipTUqKtLk0NQl1LO/e9sYO3t5cxK30Uj1w1jdSYUKvL8hi3BYwx5oK+XhORchEZbYwpE5HRQG8ryc8D5ovIbUA4ECgiTcaYk/XXjAj+fjbuvjCHORmx3P3yVi5/bC23npPJ7QuyCPL3s7o8NUCf7qvkhy9vp6rJzj2Lcvj22eN8otXSk1ixYJiIPAxUG2N+KyL3AjHGmHtOsv3NQK4x5rv97Ts3N9fk5+e7rliL1bd28Ku3dvHyphKyE8P53ZXamvF2Nc3t/OrtXby2+QiZ8WH8z9XTvWq9aFcQkU3GmNz+trOqD+a3wEIR2QcsdD5GRHJF5AmLavJKUSEBPHzlVP7v5lk0tHZy2Z/W8POVO6lv1eVRvI0xhlc2lbDg96tZubWU7543nlW3zx9x4XIqLGnBuNNIa8H0VN/awcPv7mHF+sPEhgXyo0UT+PqMFGw+1uz2RrvLGvjFm7vIK6xmZtoofn356eQkRVhdltsMtAWjATMM7ThSz3/+YwebD9cxY2w0P714oleuS+wLKhrbeOS9vbyYX0xUSAA/uDCH62aPHfGhrwEzwjkchlc3l/DgO3uoamrnK5MS+cGFOWQljty/mt6kpb2TJz87yGOrD9DR5eCmeel87/wsokIDrC7NIzRgfESTvfsXffknhbS0d/K1GSnceUEWKaN851KoJ7W0d/JsXhHLPymkurmdRZOSuHfxBNJH8N24vdGA8TE1ze089tF+nllXhMNhuHTaGG49J1NbNC7SbO/k2XXdwVLT3M7Z2fHcsSCLmWm+eWqqAeOjSuta+eunhbywoZjWji4uOC2R75w7zmffCEN1pK6VZ9Ye4vkNh2lo6+Sc7HjuuCCLGT7e56UB4+Nqm9t5Ou8QT609RF1LB1NSorh+zlgumZpMaKCuGHwyxhg2H67lyTWHeGfHUQAWTU7iW2dlaGe6kwaMArr7DF7ZVMJz64rYW95ERLA/X5+RwjWzU5mQFGl1eV6lstHO61tKeGVTCXvLm4gM9ufaOWO5cV46Y6JDrC7Pq2jAqOMYY8gvqmXFuiJWfX6U9i4HE5IiWDItmUumJPvU+JieWto7+bigklc3l/BRQSVdDsP0sdFcOTOVy6Zra68vGjCqTzXN7by1vZR/bC1lU1EtADPTRrF4chLnT0gYVmsfD0ZDWwcf7q7gnR1HWb23grYOB/ERQXxtxhiunJnC+ATtGO+PBowakOKaFlZuK+XNbaXsca5ykBEXxvkTEjh/QgIz00YRHDC8B1g6HIZdZQ18uq+Kz/ZXsvFgLe1dDhIigvjKpCQWT05idkYM/n46g+xAacCoU1Zc08JHBRV8sLuCvMJq2jsdBPrZmJISxeyMGGZnxDAzbRQRwd59M5m9s4tdpQ1sLa4jv6iWtfurqG3pHruVkxjB2dlxLJqcxPTUUSP+jlt30YBRQ9LS3sm6wmrWH6xhfWENO47U0+kwiEB6bBgTkyOZODqSScmRTEiKJDEyyJKlNupa2tlX0cS+8ib2ljeytbiOXaUNtHc5ABgdFcy8cbHMz4rjzHFxJEQGe7zGkUgDRrlUs72TLYfr2FRUy66yenaVNVBc0/rF68EBNtJiwhgbG0p6bCjJ0SHEhQc5PwKJCw8iMiRgwPOhdDkMjW0dNLR2Ut/aQXlDG2X1rZTVt1FW38aRulYKK5uPW4Y3NNCPyWOimJ4azfSx0UxLHUVSlAaKOww0YLSLXA1IWJA/Z2XFcVZW3BfP1bd2sLusgb3ljRRVt1BU3cKhqmY+2VuJvdPR634C/WwEB9gIDvAjJNAPmwhdDkOXw+Awhk6HobW9i6Y+VlfwtwmJkcGMjgrm/AnxZCVEMD4xnKyEcJKjQvSUx8towKhBiwoJYG5mLHMzY4973uEw1La0U93cTlWjnSrn58a2Tlo7umhzfrR2dOEw4Cdgswn+NsHPJoQE+BMZ4k9kcACRIQFEBvsTHxH0RavI12aFG840YJTL2WxCbHgQseFBZOtYKJ+m1+WUUm6jAaOUchsNGKWU22jAKKXcRgNGKeU2GjBKKbfRgFFKuY0GjFLKbUbcWCQRqQSKLDh0HFBlwXEHypvr8+bawLvrs6q2NGNMfH8bjbiAsYqI5A9k8JdVvLk+b64NvLs+b64N9BRJKeVGGjBKKbfRgHGd5VYX0A9vrs+bawPvrs+ba9M+GKWU+2gLRinlNhowSim30YAZJBGJEZH3RWSf83Ofa4qKSKSIHBGRR72pPhGZJiJ5IrJTRLaLyNVurmmRiBSIyH4RubeX14NE5EXn6+tFJN2d9ZxibXeJyC7nz+kDEUnzVG0Dqa/HdleIiBERr7h0rQEzePcCHxhjsoAPnI/78kvgY49U9W8Dqa8FuNEYMwlYBPy3iES7oxgR8QP+BCwGJgLXisjEEzZbCtQaY8YDfwAedEctg6xtC5BrjJkCvAI85InaTqE+RCQCuB1Y76na+qMBM3iXAk87v34auKy3jURkJpAIvOehuo7ptz5jzF5jzD7n16VABdDv3ZmDNBvYb4wpNMa0Ay84a+yr5leABeKZtVD6rc0Y85ExpsX5cB2Q4oG6Blyf0y/pDr42D9Z2Uhowg5dojCkDcH5OOHEDEbEBvwd+6OHaYAD19SQis4FA4ICb6hkDFPd4XOJ8rtdtjDGdQD0Qi/sNpLaelgL/dGtFx+u3PhGZDqQaY97yYF390km/T0JE/gUk9fLSfQPcxW3AKmNMsTv+ELugvmP7GQ08C9xkjOl9vZGh6+0HcOI9EgPZxh0GfFwRuQHIBc5xa0UnHLaX576oz/mH7A/AzZ4qaKA0YE7CGHNBX6+JSLmIjDbGlDnfoBW9bDYPmC8itwHhQKCINBljTtZf48n6EJFI4G3gfmPMOlfU1YcSILXH4xSgtI9tSkTEH4gCatxY06nUhohcQHd4n2OMsZ/4uhv1V18EMBlY7fxDlgSsFJElxhhrVyE0xujHID6Ah4F7nV/fCzzUz/Y3A496U310nxJ9ANzpgXr8gUIgw3ncbcCkE7b5f8CfnV9fA7zkoZ/VQGqbTvfpY5YFv2v91nfC9qvp7pD2aJ291mJ1AcP1g+6+gQ+Afc7PMc7nc4Enetne0wHTb33ADUAHsLXHxzQ31nQRsNf5Rr3P+dwvgCXOr4OBl4H9wAYg04M/r/5q+xdQ3uPntNLDv28nre+Ebb0mYHSogFLKbfQqklLKbTRglFJuowGjlHIbDRillNtowCil3EYDRllCRFJF5KCIxDgfj3I+9ugoZeVeGjDKEsaYYuBx4LfOp34LLDfGWLHkjHITvQ9GWUZEAoBNwJPALcB00z1aWI0QOhZJWcYY0yEiPwTeAS7UcBl59BRJWW0xUEb3YD01wmjAKMuIyDRgITAX+L5z1LcaQTRglCWcM9U9TvdI7sN0j/7+nbVVKVfTgFFWuQU4bIx53/n4MWCCiHhyIiflZnoVSSnlNtqCUUq5jQaMUsptNGCUUm6jAaOUchsNGKWU22jAKKXcRgNGKeU2/x9SW8het8dD9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"算法可视化\n",
    "可以多次运行本代码块，直观的体验分割直线存在的概率\n",
    "\"\"\"\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "# 随机取4个点为一组进行测试\n",
    "rand_s = random.randint(0,(len(xx)-4))\n",
    "rand_e = rand_s + 4\n",
    "#设置标题\n",
    "ax1.set_title('Scatter Plot')\n",
    "#设置X轴标签\n",
    "plt.xlabel('X')\n",
    "#设置Y轴标签\n",
    "plt.ylabel('Y')\n",
    "#画散点图\n",
    "ax1.scatter(xx[rand_s:rand_e],yy[rand_s:rand_e],c = 'y',marker = 'o')\n",
    "#设置图标\n",
    "plt.legend('x1')\n",
    "#画圆\n",
    "r = 0.5\n",
    "a, b = (0., 0.)\n",
    "theta = np.arange(0, 2*np.pi, 0.01)\n",
    "x = a + r * np.cos(theta)\n",
    "y = b + r * np.sin(theta)\n",
    "ax1.plot(x, y)\n",
    "#找直线\n",
    "px, py = get_exist_solve(xx[rand_s:rand_e],yy[rand_s:rand_e])\n",
    "ax1.plot([px, -px],[py, -py ])\n",
    "ax1.set_aspect(1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parent process 18322.\n",
      "Waiting for all subprocesses done...\n",
      "All subprocesses done.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"并行计算\"\"\"\n",
    "# 划分子任务\n",
    "def process(xx, yy):\n",
    "    exist_num = 0\n",
    "    for i in range(0,len(xx),4):\n",
    "        if is_exist_solve(xx[i:i+4],yy[i:i+4]):\n",
    "            exist_num += 1\n",
    "    return exist_num, i\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "# 调用的cpu数量\n",
    "cpu_num = 16\n",
    "\n",
    "# 保存结果\n",
    "result = [[0,0] for i in range(cpu_num)]\n",
    "\n",
    "# 创建进程池并启动计算\n",
    "print('Parent process %s.' % os.getpid())\n",
    "p = Pool(cpu_num)\n",
    "for i in range(cpu_num):\n",
    "    result[i] = p.apply_async(process, args=(xx[int(len(xx)/cpu_num)*i:int(len(xx)/cpu_num)*(i+1)], yy[int(len(yy)/cpu_num)*i:int(len(xx)/cpu_num)*(i+1)]))\n",
    "print('Waiting for all subprocesses done...')\n",
    "p.close()\n",
    "p.join()\n",
    "print('All subprocesses done.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "final result 0.49573707370737075\n"
     ]
    }
   ],
   "source": [
    "# 取回结果\n",
    "hit = 0\n",
    "total = 0\n",
    "for r in result:\n",
    "    hit += r.get()[0]\n",
    "    total += r.get()[1]\n",
    "    \n",
    "print('final result' ,hit/total*4)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}