github.com/cznic/mathutil@v0.0.0-20181122101859-297441e03548/example4/main.go (about)

     1  // Copyright (c) 2011 jnml. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  // +build ignore
     6  
     7  // Let QRN be the number of quadratic residues of N.  Let Q be QRN/N.  From a
     8  // sorted list of primorial products < 2^32 find "record breakers".  "Record
     9  // breaker" is N with new lowest Q.
    10  //
    11  // There are only 49 "record breakers" < 2^32.
    12  //
    13  // To run the example $ go run main.go
    14  package main
    15  
    16  import (
    17  	"fmt"
    18  	"math"
    19  	"sort"
    20  	"time"
    21  
    22  	"github.com/cznic/mathutil"
    23  	"github.com/cznic/sortutil"
    24  )
    25  
    26  func main() {
    27  	pp := mathutil.PrimorialProductsUint32(0, math.MaxUint32, 32)
    28  	sort.Sort(sortutil.Uint32Slice(pp))
    29  	var bestN, bestD uint32 = 1, 1
    30  	order, checks := 0, 0
    31  	var ixDirty uint32
    32  	m := make([]byte, math.MaxUint32>>3)
    33  	for _, n := range pp {
    34  		for i := range m[:ixDirty+1] {
    35  			m[i] = 0
    36  		}
    37  		ixDirty = 0
    38  		checks++
    39  		limit0 := mathutil.QScaleUint32(n, bestN, bestD)
    40  		if limit0 > math.MaxUint32 {
    41  			panic(0)
    42  		}
    43  		limit := uint32(limit0)
    44  		n64 := uint64(n)
    45  		hi := n64 >> 1
    46  		hits := uint32(0)
    47  		check := true
    48  		fmt.Printf("\r%10d %d/%d", n, checks, len(pp))
    49  		t0 := time.Now()
    50  		for i := uint64(0); i < hi; i++ {
    51  			sq := uint32(i * i % n64)
    52  			ix := sq >> 3
    53  			msk := byte(1 << (sq & 7))
    54  			if m[ix]&msk == 0 {
    55  				hits++
    56  				if hits >= limit {
    57  					check = false
    58  					break
    59  				}
    60  			}
    61  			m[ix] |= msk
    62  			if ix > ixDirty {
    63  				ixDirty = ix
    64  			}
    65  		}
    66  
    67  		adjPrime := ".." // Composite before
    68  		if mathutil.IsPrime(n - 1) {
    69  			adjPrime = "P." // Prime before
    70  		}
    71  		switch mathutil.IsPrime(n + 1) {
    72  		case true:
    73  			adjPrime += "P" // Prime after
    74  		case false:
    75  			adjPrime += "." // Composite after
    76  		}
    77  
    78  		if check && mathutil.QCmpUint32(hits, n, bestN, bestD) < 0 {
    79  			order++
    80  			d := time.Since(t0)
    81  			bestN, bestD = hits, n
    82  			q := float64(hits) / float64(n)
    83  			fmt.Printf(
    84  				"\r%2s #%03d %d %d %.2f %.2E %s %s %v\n",
    85  				adjPrime, order, n, hits,
    86  				1/q, q, d, time.Now().Format("15:04:05"), mathutil.FactorInt(n),
    87  			)
    88  		}
    89  	}
    90  }