{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shanks(p,g,h):      # Solves g^x=h in F_p using the Shanks algorithm.  \n",
    "    n=floor(p**0.5)+1\n",
    "    list_a=[]\n",
    "    i=0\n",
    "    temp=1\n",
    "    print 'Producing the list of powers of g...',\n",
    "    while i<=n-1:\n",
    "        list_a.append( [temp,i] )\n",
    "        temp=(temp*g)%p\n",
    "        i=i+1\n",
    "    print 'done.'                # list_a now contains [g^0,0],[g^1,1],..,[g^{n-1},n-1]\n",
    "                                 # We are keeping track of both the power AND the exponent to make our life\n",
    "                                 # easier later\n",
    "    print 'Producing the list of h*g^{-jn}...',\n",
    "    list_b=[]\n",
    "    ginv=fastpower(g,p-2,p)\n",
    "    u=fastpower(ginv,n,p)\n",
    "    i=0\n",
    "    temp=h\n",
    "    while i<=n-1:\n",
    "        list_b.append( temp)\n",
    "        temp=(temp*u)%p\n",
    "        i=i+1\n",
    "    print 'done.'               # list_b now contains hg^0,h*g^{-n},h*g^{-2n},...,h*g^{-(n-1)n}\n",
    "    print 'Starting collision algorithm for two lists of length', len(list_a)\n",
    "    \n",
    "    l=collision(list_a,list_b)  # Do the collision algorithm.  l[0]=i will be the index in the first list\n",
    "                                # and l[1]=j will be the index in the second list\n",
    "    if l=='F':\n",
    "        print \"\\n There is no solution\"\n",
    "        return\n",
    "    return l[0]+n*l[1]          # The solution is x=i+n*j.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fastpower(x,e,p):           # Compute x^e in Z/p using the fast powering algorithm\n",
    "    temp=1\n",
    "    list=[]\n",
    "    if e==0:\n",
    "        return 1\n",
    "    while e>0:\n",
    "        if e%2==0:\n",
    "            list.append(2)\n",
    "            e=e/2\n",
    "            continue\n",
    "        list.append(1)\n",
    "        e=e-1\n",
    "        continue\n",
    "    a=list.pop()\n",
    "    temp=x\n",
    "    while len(list)>0:\n",
    "        a=list.pop()\n",
    "        if a==1:\n",
    "            temp=(temp*x)%p\n",
    "            continue\n",
    "        temp=(temp*temp)%p\n",
    "    return temp\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 2, 2, 1]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "512"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fastpower(2,9,1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def collision(l,m):          #Given two lists l=[(a0,b0),(a1,b1),...] and m=[x0,x1,x2,...]\n",
    "                             #having the same length, see if some ai matches an xj.\n",
    "                             #If no, return 'F' for 'Fail'.\n",
    "    i=0                      #If yes, return the pair [i,j] giving the indices for the matching elements\n",
    "    print 'Sorting the first list...',\n",
    "\n",
    "    l.sort(key=sortFirst)  #First sort the list l so that a0<=a1<=a2<=...\n",
    "    print 'done.'               #This can take a while for very long lists.\n",
    "    print 'Testing elements of the second list...', \n",
    "    i=0                         #Now test elements of m one by one to see if they match one of the ai terms.\n",
    "    n=len(m)\n",
    "    while i<=n-1:\n",
    "        if i%100==0:\n",
    "           print '.',\n",
    "        x=m[i]\n",
    "        u=is_in(l,x)      #The function is_in decides if x matches one of the ai\n",
    "        if u=='F':\n",
    "            i=i+1\n",
    "            continue\n",
    "        print 'done'\n",
    "        return [u,i]\n",
    "    print 'done.'\n",
    "    return 'F'\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_in(l,x):          #Given l=[[a0,b0],[a1,b1],...] and x, decide if x equals one of the ai.\n",
    "                         #Use the method of \"compare to middle term\" and divide the problem in half.\n",
    "    if len(l)==0:        #If no, return 'F'.  If yes, return the bi.  \n",
    "        return 'F'       #Assumes that a0<=a1<=a2<=... from the beginning.\n",
    "    first=0           \n",
    "    last=len(l)-1        #Set up pointers \"first\" and \"last\", originally at the start and end of the list l.\n",
    "    while first<last:\n",
    "        middle=first+((last-first)//2)  #Go to roughly the middle term\n",
    "        u=l[middle][0]                  #and pull out the a-value\n",
    "        if u==x:                        #If it matches our x, return the corresponding b-value.\n",
    "            return l[middle][1]\n",
    "        if u<x:                         #If x is bigger, reset the pointers to look at the last half of the list\n",
    "            first=middle+1              #and repeat\n",
    "            continue\n",
    "        last=middle-1                   #If x is smaller, reset the pointers to look at the first half of the list.\n",
    "                                        #If we come out of the loop, we never found our x.\n",
    "    if first>last:                      #Either we exhausted all elements, or we narrowed it down to exactly one.\n",
    "        return 'F'\n",
    "    u=l[first]\n",
    "    if u[0]==x:                         #Do one last check to see if we have the x.\n",
    "        return u[1]\n",
    "    return 'F'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findprime(lower,upper,tries):\n",
    "    i=1\n",
    "    while i<=tries:\n",
    "        a=randint(lower,upper)\n",
    "        if is_prime(a):\n",
    "            return a\n",
    "        i=i+1\n",
    "    print 'None found'\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "def order(x,p):\n",
    "    i=1\n",
    "    temp=x\n",
    "    while i<=p:\n",
    "#        if i%100==0:\n",
    "#          print '\\r testing ',i,\n",
    "        if temp==1:\n",
    "            return i\n",
    "        temp=(temp*x)%p\n",
    "        i=i+1\n",
    "    print 'Error'\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35295289"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_prim_gen(p):     # Find a primitive generator for F_p\n",
    "    x=2                   # Warning: For big primes (9 or more digits) this can take a while.\n",
    "    while x<=p-1:\n",
    "        print 'Testing ',x,'...',\n",
    "        u=order(x,p)\n",
    "        print 'Order is ',u,'...',\n",
    "        if u==p-1:\n",
    "            print 'So found a generator'\n",
    "            return x\n",
    "        print 'So not a generator'\n",
    "        x=x+1\n",
    "    print 'Error'\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Testing  2 ... Order is  18198982 ... So not a generator\n",
      "Testing  3 ... Order is  72795928 ... So found a generator\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "find_prim_gen(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dlp_naive(p,g,h):   # Solve g^x=h in Z/p, given g and h.  Uses the brute force algorithm.  \n",
    "    i=0                 # i keeps track of the exponent we are trying\n",
    "    temp=1              #Start with g^0\n",
    "    while i<=p:\n",
    "        if temp==h:     #if we found it, return the exponent\n",
    "            return i\n",
    "        temp=(temp*g)%p # Otherwise multiply by one more copy of g and increase i.\n",
    "        i=i+1\n",
    "    print 'Error: no solution found'\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "83562159"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dlp_naive(a,2,55)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dlp_naive(23,5,19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ph_estimate(n):\n",
    "    i=1\n",
    "    save=0\n",
    "    while i<=n:\n",
    "        b=findprime(10^20,10^21,35)\n",
    "        if b==0:\n",
    "            continue\n",
    "        l=factor(b-1)\n",
    "        count=0\n",
    "        for a in l:\n",
    "            count=count+a[1]*a[0]\n",
    "        save=save+(b+0.0)/count\n",
    "        i=i+1\n",
    "    return (save+0.0)/n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def miller_rabin_witness(n,a):   # Determines if a is a Miller-Rabin witness for n\n",
    "    k=0                          # Returns True or False, accordingly\n",
    "    temp=n-1\n",
    "    while (temp%2)==0:           # First divide n-1 by 2 repeatedly, until we can't \n",
    "        k=k+1                    # do so anymore.\n",
    "        temp=temp/2\n",
    "    q=(n-1)/(2^k)                # Have n-1 = 2^k * q where q is odd\n",
    "    A=fastpower(a,q,n)           # A=a^q in Z/n\n",
    "    if A==1:                     # If a^q=1, we DON'T have a M-R witness\n",
    "        return False\n",
    "    i=0\n",
    "    while i<=k-1:\n",
    "        if A==n-1:               # Now test to see if A=-1 in Z/n\n",
    "            return False\n",
    "        A=(A**2)%n               # If not, square A and continue the loop\n",
    "        i=i+1\n",
    "                                 # If we get here then none of a^q,a^{2q},...,a^{2^{k-1}q}\n",
    "    return True                  # equal -1 in Z/n, so we DO have a M-R witness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_mrw(n):           # This counts the number of Miller-Rabin witnesses for n\n",
    "    i=1\n",
    "    count=0\n",
    "    while i<=n-1:\n",
    "        if miller_rabin_witness(n,i):\n",
    "            count=count+1\n",
    "        i=i+1\n",
    "    return count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "count_mrw(9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def miller_rabin_test(n,k):   # Randomly choose a number between 2 and n-1\n",
    "        i=1                   # and check if it is a M-R witness for n.\n",
    "        while i<=k:           # If yes, return it\n",
    "            a=randint(2,n-1) # If no, do it again...but stop after k tries.\n",
    "            if miller_rabin_witness(n,a):\n",
    "                print 'The number is composite: a Miller-Rabin witness is',a\n",
    "                return\n",
    "            i=i+1\n",
    "        print 'The number is probably prime'\n",
    "        return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_prime_mr(n,k):     #Test if n is prime using Miller-Rabin, k times\n",
    "    if n==2:\n",
    "        return True\n",
    "    if n==1: \n",
    "        return False\n",
    "    j=1\n",
    "    while j<=k:\n",
    "        a=randint(2,n-1)\n",
    "        if miller_rabin_witness(n,a)==True:\n",
    "            return False\n",
    "        j=j+1\n",
    "#    print \"Probably prime\"\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findaprime_mr(lower,upper,trials):\n",
    "    i=1\n",
    "    while i<=trials:\n",
    "        a=randint(lower,upper)\n",
    "        if is_prime_mr(a,10):\n",
    "            return a\n",
    "        i=i+1\n",
    "    print 'No prime found'\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pollard(N,bound):   # Perform Pollard's p-1 algorithm\n",
    "    a=2                 # to factor N, using \"bound\" as the number of steps in the loop.\n",
    "    j=2                 # See page 139 of the text.  \n",
    "    while j<=bound:\n",
    "        a=(a^j)%N\n",
    "        d=gcd(a-1,N)\n",
    "        if 1<d and d<N:\n",
    "            return d    # Have found a nontrivial factor, so return it\n",
    "        j=j+1\n",
    "    print 'No factor found'\n",
    "    return"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.7",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}