Prime sieve and factorization

[hhr]

Code: Select all

'#cmdline "-exx"
'==============================================
'The following scheme is used for the sieve.
'It has the advantage that it determines the largest and smallest divisor for each odd number.
'The smallest divisor is prime, so that only the largest divisor needs to be broken down further.

' 3*3
' 5*3-> 5*5
' 7*3-> 7*5-> 7*7
' 9*3-> 9*5-> 9*7-> 9*9
'11*3->11*5->11*7->11*9->11*11
'13*3->13*5->13*7->13*9->13*11->13*13
'15*3->15*5->15*7->15*9->15*11->15*13->15*15

'Example for 45: First 9*5 is written, then overwritten with 15*3.
'For the 15 one finds 5*3.
'For the 5 one finds 0. An unmarked odd number is prime, therefore 45=3*3*5.

'The memory requirement could be reduced by half, but I have not done this for the sake of clarity.
'==============================================

dim shared as ulong limit = 15485863
dim shared as ulong a(1 to limit),b(1 to limit)
dim as ulong i
'==============================================
declare sub Sieve
declare sub TestArray
declare sub GetPrimes(lim as ulong = limit)
declare function GetPrimeFactors(number as ulong) as string
declare sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
declare sub SavePrimeList
'==============================================

print "Please wait a few seconds..."
Sieve
cls
TestArray

for i = 3 to 45 step 2
   print i,a(i),b(i)
next i

print "Any key to continue..." : sleep

print "Prime numbers below 100:"
GetPrimes(100)

print "Any key to continue..." : sleep

for i = limit-100 to limit
   print GetPrimeFactors(i)
next i

'SavePrimeList : print "FirstMillionPrimes.txt is stored in ";exepath

'SavePrimeFactorList(2,100000) : print "PrimeFactorList.txt is stored in ";exepath

print "Any key to end..." : sleep

'==============================================

sub Sieve
   dim as ulong i,j,k
   for i = 3 to limit step 2
      for j = 3 to i step 2
         k = i * j
         if k > limit then exit for
         a(k) = i
         b(k) = j
      next
   next
end sub

sub TestArray
   dim as boolean ok = true
   dim as ulong i
   for i = 1 to limit
      if b(i) > 0 then
         if a(b(i)) > 0 then ok = false : print "b(";i;") is not prime" : sleep
      end if
   next i
   if ok = true then print "Array b() OK!"
end sub

sub GetPrimes(lim as ulong = limit)
   dim as ulong i
   print 2u
   for i = 3 to lim step 2
      if a(i) = 0 then print i
   next i
end sub

function GetPrimeFactors(number as ulong) as string
   if (number = 0) or (number = 1) then return str(number) & " = ?"
   dim as string s1,s2
   dim as ulong n = number
   dim as ubyte e2
   while (n mod 2) = 0
      n \= 2
      e2 += 1
   wend
   if e2 > 0 then s1 = "(2^" & e2 & ")"
   if n > 1 then
      while b(n) > 0
         s2 &= (b(n) & "*")
         n = a(n)
      wend
      s2 &= n
   end if
   if (len(s1) > 0) and (len(s2) > 0) then s1 &= "*"
   return str(number) & " = " & s1 & s2
end function

sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
   open "PrimeFactorList.txt" for output as #1
   dim as ulong n
   for n = n1 to n2
      print #1,GetPrimeFactors(n)
   next n
   close #1
end sub

sub SavePrimeList
   open "FirstMillionPrimes.txt" for output as #1
   dim as ulong i
   print #1,2u;
   for i = 3 to limit step 2
      if a(i) = 0 then
         print #1,
         print #1,i;
      end if
   next i
   close #1
end sub

[dodicat]

Nice.
Here is another direct method for prime factors:

Code: Select all

#cmdline "-exx"
Type primefactors
    As Ulongint f(Any)
    As Ulongint number
    Declare Operator Cast() As String
End Type

Operator primefactors.cast() As String'custom print
Dim As String s=Str(number)+"  =  "
Dim As String op
For m As Long=1 To Ubound(f)
    If m<Ubound(f) Then op="*" Else op=""
    s+=Str(f(m))+op
Next m
Dim As String d
If Ubound(f)=1 Then d=String(5,"-")+">(PRIME)" Else d=""
Return s + d
End Operator

Function getprimefactors(number As Ulongint,p As primefactors) As Integer
    Dim As Integer i
    p.number=number
    Redim p.f(1 to 64)
    #macro fill(x)
    i+=1
    p.f(i)=x
    #endmacro
    #macro Eliminate(x)
    While number Mod x=0
        number=number\x
        fill(x)
        check=check*x
        limit=Sqr(number)+1
    Wend
    #endmacro
    Dim As Ulongint original=number,divisor=0ull,limit=Sqr(number)+1,check=1ull
    
    Eliminate(2)
    Eliminate(3)
    
    While divisor<limit
        divisor=divisor+6 
        Eliminate((divisor-1))
        Eliminate((divisor+1))
    Wend
    If number>1 Then 
        fill(number)
        check=check*number
    End If
    number=original
    redim preserve p.f(1 to i)
    Return check=original
End Function

Dim As primefactors pf
Print "Prime numbers under 100"
For n As Long=2 To 100
    If Not GetPrimeFactors(n,pf) Then Print "Error" Else If Ubound(pf.f)=1 Then Print pf 
Next n
Print "press a key"
Sleep
Print
Dim  As Ulong limit = 15485863
For i As Ulongint= limit-100 To limit
    Color 15
    Var gpf= GetPrimeFactors(i,pf)
    If Not gpf Then Print "Error"
    If Instr(pf,"PRIME") Then Color 2 
    Print pf
Next i


Sleep 

[hhr]

I have read that one can check whether the number of prime factors is even or odd.
This should be pseudorandom.
In your example, I only needed to change line 53 for a short test, thank you.

Code: Select all

'#cmdline "-exx"
Type primefactors
    As Ulongint f(Any)
    As Ulongint number
    Declare Operator Cast() As String
End Type

Operator primefactors.cast() As String'custom print
Dim As String s=Str(number)+"  =  "
Dim As String op
For m As Long=1 To Ubound(f)
    If m<Ubound(f) Then op="*" Else op=""
    s+=Str(f(m))+op
Next m
Dim As String d
If Ubound(f)=1 Then d=String(5,"-")+">(PRIME)" Else d=""
Return s + d
End Operator

Function getprimefactors(number As Ulongint,p As primefactors) As Integer
    Dim As Integer i
    p.number=number
    Redim p.f(1 to 64)
    #macro fill(x)
    i+=1
    p.f(i)=x
    #endmacro
    #macro Eliminate(x)
    While number Mod x=0
        number=number\x
        fill(x)
        check=check*x
        limit=Sqr(number)+1
    Wend
    #endmacro
    Dim As Ulongint original=number,divisor=0ull,limit=Sqr(number)+1,check=1ull
    
    Eliminate(2)
    Eliminate(3)
    
    While divisor<limit
        divisor=divisor+6 
        Eliminate((divisor-1))
        Eliminate((divisor+1))
    Wend
    If number>1 Then 
        fill(number)
        check=check*number
    End If
    number=original
    redim preserve p.f(1 to i)
    Return i ' !!! Changed !!!
End Function

Dim Shared As primefactors pf

Function PrimeFactorPrng As Ulongint
   Dim as Ubyte i
   Dim As Ulongint b
   Static As Ulongint n = -1 ' Seed
   For i = 0 To 63
      n += 1
      if n < 2 then n = 2
      If (getprimefactors(n,pf) Mod 2) = 1 Then b = Bitset(b,i)
   Next i
   Return b
End Function

Do
   Print Bin(PrimeFactorPrng,64);
Loop Until Len(Inkey)

Print : Print

'Chdir "Path to PractRand's RNG_test.exe" '' PractRand Version 0.94
Dim As String cs = "RNG_test stdin64"
cs &= " -tlmin 1KB"
Open Pipe cs For Binary Access Write As #1
Do
   Put #1,,PrimeFactorPrng
Loop

[hhr]

Referring to the first example:

Code: Select all

'#cmdline "-exx"
'==============================================
'Because only odd numbers have to be processed,
'the memory requirement of the arrays can be halved with (limit-1)\2.
'Now only the values for the odd numbers are stored
'and the indices must be adjusted for all array accesses.
'==============================================
dim shared as ulong limit = 15485863
dim shared as ulong a(1 to ((limit-1) shr 1)), b(1 to ((limit-1) shr 1)) ' (limit-1) shr 1 = (limit-1)\2
dim as ulong i
'==============================================
declare sub Sieve
declare sub TestArray
declare sub GetPrimes(lim as ulong = limit)
declare function GetPrimeFactors(number as ulong) as string
declare sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
declare sub SavePrimeList
'==============================================

print "Please wait a few seconds..."
Sieve
cls
TestArray

for i = 3 to 45 step 2
   print i, a((i-1) shr 1), b((i-1) shr 1)
next i

print "Any key to continue..." : sleep

print "Prime numbers below 100:"
GetPrimes(100)

print "Any key to continue..." : sleep

for i = limit-100 to limit
   print GetPrimeFactors(i)
next i

'SavePrimeList : print "FirstMillionPrimes2.txt is stored in ";exepath

'SavePrimeFactorList(2,100000) : print "PrimeFactorList2.txt is stored in ";exepath

print "Any key to end..." : sleep

'==============================================

sub Sieve
   dim as ulong i,j,k,index
   for i = 3 to limit\3 step 2
      for j = 3 to i step 2
         k = i * j
         if k > limit then exit for
         index = (k-1) shr 1
         a(index) = i
         b(index) = j
      next
   next
end sub

sub TestArray
   dim as boolean ok = true
   dim as ulong i,index
   for i = 3 to limit step 2
      index = (i-1) shr 1
      if b(index) > 0 then
         if a((b(index)-1) shr 1) > 0 then ok = false : print "b(";i;") is not prime" : sleep
      end if
   next i
   if ok = true then print "Array b() OK!"
end sub

sub GetPrimes(lim as ulong = limit)
   dim as ulong i
   print 2u
   for i = 3 to lim step 2
      if a((i-1) shr 1) = 0 then print i
   next i
end sub

function GetPrimeFactors(number as ulong) as string
   if (number = 0) or (number = 1) then return str(number) & " = ?"
   dim as string s1,s2
   dim as ulong n = number
   dim as ubyte e2
   while (n mod 2) = 0
      n \= 2
      e2 += 1
   wend
   if e2 > 0 then s1 = "(2^" & e2 & ")"
   if n > 1 then
      while b((n-1) shr 1) > 0
         s2 &= (str(b((n-1) shr 1)) & "*")
         n = a((n-1) shr 1)
      wend
      s2 &= str(n)
   end if
   if (len(s1) > 0) and (len(s2) > 0) then s1 &= "*"
   return str(number) & " = " & s1 & s2
end function

sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
   open "PrimeFactorList2.txt" for output as #1
   dim as ulong n
   for n = n1 to n2
      print #1,GetPrimeFactors(n)
   next n
   close #1
end sub

sub SavePrimeList
   open "FirstMillionPrimes2.txt" for output as #1
   dim as ulong i
   print #1,2u;
   for i = 3 to limit step 2
      if a((i-1) shr 1) = 0 then
         print #1,
         print #1,i;
      end if
   next i
   close #1
end sub

[hhr]

I have transferred the example from this page to FreeBASIC:

Code: Select all

dim as ulong limit = 50
dim as ulong s(2 to limit)
dim as ulong i, p

for p = 2 to limit
   if p * p > limit then exit for
   if s(p) = 0 then
      for i = p * p to limit step p
         s(i) = 1
      next
   end if
next

for p = 2 to limit
   if s(p) = 0 then print p
next

sleep

This resulted in the following example.
The sieve only needs a quarter of the array size of the first example.
And because it overwrites array entries less frequently, it is about twice as fast as this one.
However, the output is not sorted:

Code: Select all

'#cmdline "-exx"
dim shared as ulong limit = 15485863
dim shared as ulong a(1 to (limit-1)\2)
dim as ulong i
'==============================================
declare sub SieveOfEratosthenes
declare sub GetPrimes(lim as ulong = limit)
declare function GetPrimeFactors(number as ulong) as string
declare sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
declare sub SavePrimeList
'==============================================

print "Please wait a few seconds..."
SieveOfEratosthenes

for i = 3 to 45 step 2
   print i,a(i)
next i

print "Any key to continue..." : sleep

print "Prime numbers below 100:"
GetPrimes(100)

print "Any key to continue..." : sleep

for i = limit-100 to limit
   print GetPrimeFactors(i)
next i

'SavePrimeList : print "FirstMillionPrimes3.txt is stored in ";exepath

'SavePrimeFactorList(2,100000) : print "PrimeFactorList3.txt is stored in ";exepath

print "Any key to end..." : sleep

'=================================================

sub SieveOfEratosthenes
   dim as ulong i, j
   for i = 3 to limit step 2
      if i * i > limit then exit for
      if a((i-1)\2) = 0 then
         for j = i * i to limit step 2 * i
            a((j-1)\2) = i ' i instead of 1
         next
      end if
   next
end sub

sub GetPrimes(lim as ulong = limit)
   dim as ulong i
   print 2u
   for i = 3 to lim step 2
      if a((i-1)\2) = 0 then print i
   next i
end sub

function GetPrimeFactors(number as ulong) as string
   if (number = 0) or (number = 1) then return str(number) & " = ?"
   dim as string s1,s2
   dim as ulong n = number
   dim as ubyte e2
   while (n mod 2) = 0
      n \= 2
      e2 += 1
   wend
   if e2 > 0 then s1 = "(2^" & e2 & ")"
   if n > 1 then
      while a((n-1)\2) > 0
         s2 &= (str(a((n-1)\2)) & "*") ' s2 = str(a((n-1)\2)) & "*" & s2
         n \= a((n-1)\2)
      wend
      s2 &= str(n)
   end if
   if (len(s1) > 0) and (len(s2) > 0) then s1 &= "*"
   return str(number) & " = " & s1 & s2
end function

sub SavePrimeFactorList(n1 as ulong, n2 as ulong)
   open "PrimeFactorList3.txt" for output as #1
   dim as ulong n
   for n = n1 to n2
      print #1,GetPrimeFactors(n)
   next n
   close #1
end sub

sub SavePrimeList
   open "FirstMillionPrimes3.txt" for output as #1
   dim as ulong i
   print #1,2u;
   for i = 3 to limit step 2
      if a((i-1)\2) = 0 then
         print #1,
         print #1,i;
      end if
   next i
   close #1
end sub

Edit:
After a small change, each array element is written at most once.
Because only the smallest prime factors are now stored, this results in a sorted output:

Code: Select all

sub SieveOfEratosthenes
   dim as ulong i, j
   for i = 3 to limit step 2
      if i * i > limit then exit for
      if a((i-1)\2) = 0 then
         for j = i * i to limit step 2 * i
            if a((j-1)\2) = 0 then a((j-1)\2) = i ' This line changed
         next
      end if
   next
end sub

[srvaldez]

hi hhr
now all you need is to implement the Quadratic sieve https://en.wikipedia.org/wiki/Quadratic_sieve
that way it would be feasible to factor integers up to 100 digits