Cowcave/website/source/cyclop.cpp at main · JonathanCrombie/Cowcave

1150 lines (906 loc) · 31 KB
 * This source distribution is placed into the public domain by its author, Jonathan Crombie.
 * You may use it for whatever purpose, without charge, and without having to notify anyone.
 * I disclaim any responsibility for any errors or lack of fitness for any purpose.
 * cyclop Version 1.4.4  Oct. 13, 2022
 * (Oct 13, 2022. Smarter ComputeLM. Eliminate mpz_pow() call from within for loops.)
 * (May 24, 2018. Added capability for ComputeAllFactors() to add 1 and N as factors.)
 * (Jun 30, 2016. Removed the pari dependency.)
 * (Sep 22, 2014. Improved EulerTotient() -- skip computing all composite factors.)
 * (Changes from version 1.03:
 *  -- Allow perfect powers for base.)
 * (Changes from version 1.02:
 *  -- Some optimizations.  Instead of trial factoring up to sqrt(n) for RemoveSquares(),
 *     ComputeAllFactors() etc., now only compute the prime factorization and then permute
 *     factors to compute the remaining composite factors.)
 * (Changes from Version 1.01:
 *  -- minor tweak, change pari_init() to not initialize a large number of prime differences.
 *     this speeds up initialization quite a bit.)
 * (Changes from Version 1.00:
 *  -- fix ComputePrimitive() not properly implemented.
 *  -- fix algorithm fail if exponent == 1
 *  -- don't allow exponents < 1 )
 * To build try:
 * g++ cyclop.cpp -lgmp -o cyclop
 * (You will need to have g++ and the gmp library already installed)
 * (For linux bash shell for Windows 10, try:
    sudo apt-get install g++ libgmp3-dev                           )
 * This source file is currently located at:
 * http://myfactors.mooo.com/source/cyclop.cpp
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <values.h>
#include <gmp.h>
#include <string.h>
#include <ctype.h>
struct factor_info {
mpz_t           the_factor;
long            occurrences;
long            digits;
char            factor_status;
struct factor_infos {
long                 count;
struct factor_info*  the_factors;
void ConvertPerfectPowers( mpz_t, mpz_t, mpz_t, mpz_t, mpz_t );
void ComputeLMPrimitives( mpz_t, mpz_t, mpz_t, char, mpz_t, mpz_t );
void ComputePrimitive( mpz_t, mpz_t, char, mpz_t );
void ComputeLongForm( mpz_t, mpz_t, char, mpz_t );
int IsLM( mpz_t, mpz_t, char );
int ComputeLM( mpz_t, mpz_t, mpz_t, char, mpz_t, mpz_t );
int ComputeLucasCD( mpz_t, mpz_t**, long*, mpz_t**, long* );
char quickprimecheck( mpz_t );
void TFDivideOut( mpz_t, long, char*, mpz_t, struct factor_infos* );
void Factorize( mpz_t, struct factor_infos* );
void PermuteFactors( struct factor_infos*, struct factor_infos*, long, mpz_t );
void ComputeAllFactors( mpz_t, struct factor_infos*, int = 0 );
long ComputeOccurrences( mpz_t, mpz_t );
int Mobius( mpz_t );
void EulerTotient( mpz_t, mpz_t );
void Init_Factor_Infos( struct factor_infos* );
void AddFactorInfo( struct factor_infos*, mpz_t, long, char, long );
void Cleanup_Factor_Infos( struct factor_infos* );
long NumberOfDigits( mpz_t );
void mpz_pow( mpz_t, mpz_t, mpz_t );
int factor_info_cmpfunc( const void*, const void* );
unsigned char wheel7[48] = { 2, 4, 2, 4, 6, 2, 6, 4, 2, 4,
int main(int argc, char* argv[] )
  if ( argc != 4 ) {
    printf("\n");
    printf("cyclop Version 1.4.4\n\n");
    printf("Usage: \"cyclop base exponent +|-\"");
    printf("\n\n");
    return -1;
  mpz_t base;
  mpz_init( base );
  mpz_set_str( base, argv[1], 10 );
  if ( mpz_cmp_ui( base, 2 ) < 0 ) {
    printf( "\n\nbase must be >= 2.  Aborting.\n\n" );
    mpz_clear( base );
    return -1;
  mpz_t exponent;
  mpz_init( exponent );
  mpz_set_str( exponent, argv[2], 10 );
  if ( mpz_cmp_ui( exponent, 1 ) < 0 ) {
    printf( "\n\nexponent cannot be less than 1.  Aborting.\n\n" );
    mpz_clear( exponent );
    mpz_clear( base );
    return -1;
  char c0 = argv[3][0];
  if ( c0 != '+' && c0 != '-' ) {
    printf("\n\nThird argument must be either \'+\' or \'-\'.\n\n");
    mpz_clear( exponent );
    mpz_clear( base );
    return -1;
  mpz_t ConvertedBase;
  mpz_init( ConvertedBase );
  mpz_t ConvertedExponent;
  mpz_init( ConvertedExponent );
  mpz_t SqrFreeBase;
  mpz_init( SqrFreeBase );
  ConvertPerfectPowers( base, exponent, ConvertedBase, ConvertedExponent, SqrFreeBase );
  if ( IsLM( SqrFreeBase, ConvertedExponent, c0 ) ) {
    mpz_t LPrimitive;
    mpz_init( LPrimitive );
    mpz_t MPrimitive;
    mpz_init( MPrimitive );
    ComputeLMPrimitives( ConvertedBase, SqrFreeBase, ConvertedExponent, c0, LPrimitive, MPrimitive );
    gmp_printf( "L-Primitive=%Zd\n", LPrimitive );
    gmp_printf( "M-Primitive=%Zd\n", MPrimitive );
    mpz_clear( MPrimitive );
    mpz_clear( LPrimitive );
    mpz_t Primitive;
    mpz_init( Primitive );
    ComputePrimitive( ConvertedBase, ConvertedExponent, c0, Primitive );
    gmp_printf( "Primitive=%Zd\n", Primitive );
    mpz_clear( Primitive );
  mpz_clear( ConvertedBase );
  mpz_clear( ConvertedExponent );
  mpz_clear( SqrFreeBase );
  mpz_clear( exponent );
  mpz_clear( base );
  return 0;
// Need to fix base, exponent if base is a perfect power
// Compute the square-free part of base as well, since it's really easy to do.
void ConvertPerfectPowers( mpz_t base, mpz_t exponent, mpz_t ConvertedBase, mpz_t ConvertedExponent, mpz_t SqrFreeBase ) {
  struct factor_infos PrimeFactorization;
  Init_Factor_Infos( &PrimeFactorization );
  Factorize( base, &PrimeFactorization );
  mpz_t gcd;
  mpz_init_set_ui( gcd, PrimeFactorization.the_factors[0].occurrences );
  for ( i = 1; i < PrimeFactorization.count && mpz_cmp_ui( gcd, 1 ) != 0; i++ )
    mpz_gcd_ui( gcd, gcd, PrimeFactorization.the_factors[i].occurrences );
  long power = mpz_get_d( gcd );
  mpz_set_ui( ConvertedBase, 1 );
  mpz_mul_ui( ConvertedExponent, exponent, power );
  mpz_t tempz1;
  mpz_init( tempz1 );
  mpz_set_ui( SqrFreeBase, 1 );
  for ( i = 0; i < PrimeFactorization.count; i++ ) {
    long reducedpower = PrimeFactorization.the_factors[i].occurrences / power;
    mpz_pow_ui( tempz1, PrimeFactorization.the_factors[i].the_factor, reducedpower );
    mpz_mul( ConvertedBase, ConvertedBase, tempz1 );
    if ( reducedpower % 2 )
      mpz_mul( SqrFreeBase, SqrFreeBase, PrimeFactorization.the_factors[i].the_factor );
  mpz_clear( tempz1 );
  mpz_clear( gcd );
  Cleanup_Factor_Infos( &PrimeFactorization );
// Algorithm from "Factorizations of b^n ± 1, b = 2, 3, 5, 6, 7, 10, 11, 12 Up to High Powers", Third Edition, page lxxi,
//  John Brillhart, D. H. Lehmer, J. L. Selfridge, Bryant Tuckerman, and S. S. Wagstaff, Jr.
// Correction information for bases with odd square factors provided by user "jyb" on:
//  http://www.mersenneforum.org/showpost.php?p=186049&postcount=5
// Assumption: base is not a perfect power
void ComputeLMPrimitives( mpz_t base, mpz_t sqrfreebase, mpz_t exponent, char c0, mpz_t LPrimitive, mpz_t MPrimitive ) {
  mpz_t m;  // odd part of exponent
  mpz_init_set( m, exponent );
  while ( mpz_even_p( m ) )
    mpz_div_2exp( m, m, 1 );
  struct factor_infos Factor_Infos;
  Init_Factor_Infos( &Factor_Infos );
  ComputeAllFactors( m, &Factor_Infos, 1 );
  mpz_t LPrimNumerator;
  mpz_init_set_ui( LPrimNumerator, 1 );
  mpz_t LPrimDenominator;
  mpz_init_set_ui( LPrimDenominator, 1 );
  mpz_t MPrimNumerator;
  mpz_init_set_ui( MPrimNumerator, 1 );
  mpz_t MPrimDenominator;
  mpz_init_set_ui( MPrimDenominator, 1 );
  mpz_t gcd;
  mpz_init( gcd );
  mpz_init( d );
  long i = 0;
  for ( ; i < Factor_Infos.count; i++ ) {
    mpz_set( d, Factor_Infos.the_factors[i].the_factor );
    mpz_gcd( gcd, sqrfreebase, d );
    if ( mpz_cmp_ui( gcd, 1 ) != 0 )
      continue;
    mpz_t L;
    mpz_init( L );
    mpz_t M;
    mpz_init( M );
    mpz_t cofactor;
    mpz_init( cofactor );
    mpz_divexact( cofactor, exponent, d );
    ComputeLM( base, sqrfreebase, cofactor, c0, L, M );
    int jacobi = mpz_jacobi( sqrfreebase, d );
    int mobius = Mobius( d );
    if ( mobius != 0 ) {
      mpz_mul( mobius == 1 ? LPrimNumerator : LPrimDenominator, mobius == 1 ? LPrimNumerator : LPrimDenominator, jacobi == 1 ? L : M );
      mpz_mul( mobius == 1 ? MPrimNumerator : MPrimDenominator, mobius == 1 ? MPrimNumerator : MPrimDenominator, jacobi == 1 ? M : L );
    mpz_clear( cofactor );
    mpz_clear( M );
    mpz_clear( L );
  mpz_divexact( LPrimitive, LPrimNumerator, LPrimDenominator );
  mpz_divexact( MPrimitive, MPrimNumerator, MPrimDenominator );
  mpz_clear( d );
  mpz_clear( gcd );
  mpz_clear( MPrimDenominator );
  mpz_clear( MPrimNumerator );
  mpz_clear( LPrimDenominator );
  mpz_clear( LPrimNumerator );
  Cleanup_Factor_Infos( &Factor_Infos );
  mpz_clear( m );
// Algorithm from A. Kruppa's phi.c found at <http://www.asahi-net.or.jp/~KC2H-MSM/mathland/matha1/matha152x.htm>
void ComputePrimitive( mpz_t base, mpz_t exponent, char c0, mpz_t Primitive ) {
  if ( mpz_cmp_ui( exponent, 1 ) == 0 ) {
    ComputeLongForm( base, exponent, c0, Primitive );
    return;
  mpz_t phi_exponent;
  mpz_init_set( phi_exponent, exponent );
  if ( c0 == '+' )
    mpz_mul_2exp( phi_exponent, phi_exponent, 1 ); // multiply by 2^1
  mpz_t Rp;
  mpz_init_set_ui( Rp, 1 );
  mpz_t Rq;
  mpz_init_set_ui( Rq, 1 );
  mpz_t tempz1;
  mpz_init( tempz1 );
  struct factor_infos Factor_Infos;
  Init_Factor_Infos( &Factor_Infos );
  ComputeAllFactors( phi_exponent, &Factor_Infos, 1 );
  long i = 0;
  for ( ; i < Factor_Infos.count; i++ ) {
    ComputeLongForm( base, Factor_Infos.the_factors[i].the_factor, '-', tempz1 );
    mpz_t cofactor;
    mpz_init( cofactor );
    mpz_divexact( cofactor, phi_exponent, Factor_Infos.the_factors[i].the_factor );
    switch ( Mobius( cofactor ) ) {
      case 1:
        mpz_mul( Rp, Rp, tempz1 );
        break;
      case -1:
        mpz_mul( Rq, Rq, tempz1 );
        break;
      case 0:
      default:
        break;
    mpz_clear( cofactor );
  mpz_divexact( Primitive, Rp, Rq );
  Cleanup_Factor_Infos( &Factor_Infos );
  mpz_clear( tempz1 );
  mpz_clear( Rq );
  mpz_clear( Rp );
  mpz_clear( phi_exponent );
// Compute the value base^exponent + c0
void ComputeLongForm( mpz_t base, mpz_t exponent, char c0, mpz_t the_number ) {
  if ( mpz_cmp_ui( base, 2 ) == 0 ) {
    mpz_t number1;
    mpz_init_set_ui( number1, 1 );
    mp_bitcnt_t bitcount = mpz_get_ui( exponent );
    mpz_mul_2exp( the_number, number1, bitcount );
    if ( c0 == '-' )
      mpz_sub( the_number, the_number, number1 );
      mpz_add( the_number, the_number, number1 );
    mpz_clear( number1 );
    return;
  if ( mpz_cmp_ui( exponent, ULONG_MAX ) > 0 ) {
    printf("Error: exponent must be <= %lu\n", ULONG_MAX );
    mpz_t number1;
    mpz_init_set_ui( number1, 1 );
    unsigned long ulexp = mpz_get_ui( exponent ); 
    mpz_pow_ui( the_number, base, ulexp );
    if ( c0 == '-' )
      mpz_sub( the_number, the_number, number1 );
      mpz_add( the_number, the_number, number1 );
    mpz_clear( number1 );
// Returns 1 if base^exponent + c0 has an Aurifeuillian factorization else 0.
// Assumption: base is not a perfect power and is square-free
int IsLM( mpz_t base, mpz_t exponent, char c0 ) {
  // Not handled values
  if ( mpz_cmp_ui( base, 2 ) < 0 ||
       mpz_cmp_ui( exponent, 2 ) < 0 ||
       (c0 != '-' && c0 != '+') ) {
    return 0;
  {  // check that we have the correct c0.  c0 should be -1 if base == 1 mod 4, else c0 should be 1
    mpz_t r;
    mpz_init( r );
    mpz_mod_ui( r, base, 4 );
    if (  ( (c0 == '+') && mpz_cmp_ui( r, 1 ) == 0) ||
          ( (c0 == '-') && mpz_cmp_ui( r, 1 ) != 0) ) {
      mpz_clear( r );
      return 0;
    mpz_clear( r );
  if ( !mpz_divisible_p( exponent, base ) )  // not divisible.  L,M do not exist
    return 0;
  mpz_init( h );
  mpz_divexact( h, exponent, base );
  if ( mpz_even_p( h ) ) { // h is even.  L,M do not exist
    mpz_clear( h );
    return 0;
  mpz_clear( h );
  return 1;
// Compute C and D polynomial coefficients.  Code implemented based on
// Richard Brent's "On computing factors of cyclotomic polynomials" Algorithm L.
int ComputeLucasCD( mpz_t n, mpz_t** C, long* C_count, mpz_t** D, long* D_count ) {
  mpz_t tempz1;
  mpz_init( tempz1 );
  mpz_t tempz2;
  mpz_init( tempz2 );
  mpz_t nmod4;
  mpz_init( nmod4 );
  mpz_mod_ui( nmod4, n, 4 );
  mpz_t nmod8;
  mpz_init( nmod8 );
  mpz_mod_ui( nmod8, n, 8 );
  mpz_t nprime;
  mpz_init( nprime );
  if ( mpz_cmp_ui( nmod4, 1 ) == 0 )
    mpz_set( nprime, n );
    mpz_mul_ui( nprime, n, 2 );
  mpz_init( s );
  if ( mpz_cmp_ui( nmod4, 3 ) == 0 )
    mpz_set_si( s, -1 );
    mpz_set_ui( s, 1 );
  mpz_t sprime;
  mpz_init( sprime );
  if ( mpz_cmp_ui( nmod8, 5 ) == 0 )
    mpz_set_si( sprime, -1 );
    mpz_set_ui( sprime, 1 );
  mpz_init( d );
  EulerTotient(nprime, tempz1 );
  mpz_divexact_ui( d, tempz1, 2 );
  // Step 1
  long ld = mpz_get_ui( d );
  mpz_t* q = (mpz_t*) calloc( ld + 1, sizeof( mpz_t ) );
  for ( k = 0; k <= ld; k++ )
    mpz_init_set_ui( q[k], 0 );
  mpz_t mpz_k;
  mpz_init( mpz_k );
  for ( k = 1; k <= ld; k++ ) {
    mpz_set_ui( mpz_k, k );
    if ( k % 2 )
      mpz_set_si( q[k], mpz_jacobi( n, mpz_k ) );
      mpz_t gprime_k;
      mpz_init( gprime_k );
      mpz_gcd( gprime_k, mpz_k, nprime );
      mpz_t tempz1;
      mpz_init( tempz1 );
      mpz_divexact( tempz1, nprime, gprime_k );
      int mobius_value = Mobius( tempz1 );
      mpz_t mpz_totient_value;
      mpz_init( mpz_totient_value );
      EulerTotient( gprime_k, mpz_totient_value );
      int cosine_value = 0;
        mpz_sub_ui( tempz1, n, 1 );
        mpz_mul( tempz1, tempz1, mpz_k );
        mpz_t mod8;
        mpz_init( mod8 );
        mpz_mod_ui( mod8, tempz1, 8 );
        // we won't check for odd values since k is even and therefore mod8 is also even
        if ( mpz_cmp_ui( mod8, 0 ) == 0 )
          cosine_value = 1;
        else if ( mpz_cmp_ui( mod8, 2 ) == 0 )
          cosine_value = 0;
        else if ( mpz_cmp_ui( mod8, 4 ) == 0 )
          cosine_value = -1;
        else // mod8 == 6
          cosine_value = 0;
        mpz_clear( mod8 );        
      mpz_mul_si( tempz1, mpz_totient_value, mobius_value );
      mpz_mul_si( tempz1, tempz1, cosine_value );
      mpz_set( q[k], tempz1 );
      mpz_clear( mpz_totient_value );
      mpz_clear( tempz1 ); 
      mpz_clear( gprime_k );
  mpz_t* gamma = (mpz_t*) calloc( ld + 1, sizeof( mpz_t ) );
  for ( k = 0; k <= ld; k++ )
    mpz_init_set_ui( gamma[k], 0 );
  mpz_t* delta = (mpz_t*) calloc( ld, sizeof( mpz_t ) );
  for ( k = 0; k < ld; k++ )
    mpz_init_set_ui( delta[k], 0 );
  // Step 2
  mpz_set_ui( gamma[0], 1 );
  mpz_set_ui( delta[0], 1 );
  // Step 3
  int d_is_odd = ld % 2;
  long max_gamma = d_is_odd ? ( ld - 1 ) / 2 : ld / 2;
  long max_delta = d_is_odd ? ( ld - 1 ) / 2 : ( ld - 2 ) / 2;
  for ( k = 1; k <= max_gamma; k++ ) {
    // compute gamma[k]
    long j = 0;
    for ( j = 0; j <= k - 1; j++ ) {
      long q_index = 2 * k - 2 * j - 1;
      mpz_mul( tempz1, n, q[q_index] );
      mpz_mul( tempz1, tempz1, delta[j] );
      mpz_add( gamma[k], gamma[k], tempz1 );
      q_index++;
      mpz_mul( tempz1, q[q_index], gamma[j] );
      mpz_sub( gamma[k], gamma[k], tempz1 );
    mpz_divexact_ui( gamma[k], gamma[k], 2 * k );
    if ( k > max_delta )
      continue;
    // compute delta[k]
    for ( j = 0; j <= k - 1; j++ ) {
      long q_index = 2 * k + 1 - 2 * j;
      mpz_mul( tempz1, q[q_index], gamma[j] );
      mpz_add( delta[k], delta[k], tempz1 );
      q_index--;
      mpz_mul( tempz1, q[q_index], delta[j] );
      mpz_sub( delta[k], delta[k], tempz1 );
    mpz_add( delta[k], delta[k], gamma[k] );
    mpz_divexact_ui( delta[k], delta[k], 2 * k + 1 );
  // Step 4
  for ( k = max_gamma + 1; k <= ld; k++ )
    mpz_set( gamma[k], gamma[ld-k] );
  // Step 5
  for ( k = max_delta + 1; k < ld; k++ )
    mpz_set( delta[k], delta[ld-1-k] );
  // Pass back
  *C = gamma;
  *C_count = ld + 1;
  *D = delta;
  *D_count = ld;
  // clean up
  mpz_clear( mpz_k );
  if ( q != NULL ) {
    for ( k = 0; k <= ld; k++ )
      mpz_clear( q[k] );
    free( q );
    q = NULL;
  mpz_clear( d );
  mpz_clear( sprime );
  mpz_clear( s );
  mpz_clear( nprime );
  mpz_clear( nmod8 );
  mpz_clear( nmod4 );
  mpz_clear( tempz2 );
  mpz_clear( tempz1 );
  return 0;
// Computes the Aurifeuillian factors L,M for a^x +/- 1.  Returns 0 on success.
// Assumption: base is not a perfect power
int ComputeLM( mpz_t base, mpz_t sqrfreebase, mpz_t exponent, char c0, mpz_t L, mpz_t M ) {
  // Not handled values
  if ( mpz_cmp_ui( base, 2 ) < 0 ||
       mpz_cmp_ui( exponent, 2 ) < 0 ||
       (c0 != '-' && c0 != '+') ) {
    return -1;
  {  // check that we have the correct c0.  c0 should be -1 if sqrfreebase == 1 mod 4, else c0 should be 1
    mpz_t r;
    mpz_init( r );
    mpz_mod_ui( r, sqrfreebase, 4 );
    if ( (c0 == '+' && mpz_cmp_ui( r, 1 ) == 0) ||
         (c0 == '-' && mpz_cmp_ui( r, 1 ) != 0) ) {
      mpz_clear( r );
      return -1;
    mpz_clear( r );
  if ( !mpz_divisible_p( exponent, sqrfreebase ) )  // not divisible.  L,M do not exist
    return -1;
  unsigned long h = 0;
    mpz_t temph;
    mpz_init( temph );
    mpz_divexact( temph, exponent, sqrfreebase );
    h = mpz_get_ui( temph );
    mpz_clear( temph );
    if ( (h % 2) == 0 ) // h is even.  L,M do not exist
      return -1;
  mpz_t* C = NULL;
  long C_count = 0;
  mpz_t* D = NULL;
  long D_count = 0;
  if ( ComputeLucasCD( sqrfreebase, &C, &C_count, &D, &D_count ) != 0 ) {
    gmp_printf("Error: Could not compute Lucas C,D polynomial coefficients for base (%Zd). Aborting.\n", sqrfreebase );
    return -1;
  mpz_t tempz1;
  mpz_init( tempz1 );
  unsigned long k = (h + 1) / 2;
  // Both Polynomial A & B are initially computed as follows:
  //   coef[0] * base^(h*0) + coef[1] * base^(h*1) + coef[2] * base^(h*2) + ...
  // = coef[0] * (base^h)^0 + coef[1] * (base^h)^1 + coef[2] * (base^h)^2 + ...
  // = coef[0] * 1 + coef[1] * base^h + coef[2] * base^h * base^h + ...
  // B polynomial is then multiplied by additional terms
  //   L = A - B
  //   M = A + B
  // NOTE: Reading coefficients from beginning to end is exactly the same as reading from the end to the beginning.
  // eg. for base 31, C = { 1,15,43,83,125,151,169,173,173,169,151,125,83,43,15,1 } and D = { 1,5,11,19,25,29,31,31,31,29,25,19,11,5,1 }
  // Also, D_count = C_count - 1
 // Compute Polynomial A with the Lucas 'C' coefficients
  mpz_init( A );
    mpz_t base_h; // base^h
    mpz_init( base_h );
    mpz_pow_ui( base_h, base, h );
    mpz_t base_raised; // base_h^i
    mpz_init( base_raised );
    for ( i = 0; i < C_count; i++ ) { // each iteration computes one term and then adds to A
      if ( i == 0 )
        mpz_set_ui( base_raised, 1 );
      else if ( i == 1 )
        mpz_set( base_raised, base_h );
        mpz_mul( base_raised, base_raised, base_h );
      mpz_mul( tempz1, base_raised, C[i] );
      mpz_add( A, A, tempz1 );
    mpz_clear( base_raised );
    mpz_clear( base_h );
  // Compute Polynomial B with Lucas 'D' coefficients
  mpz_init( B );
    mpz_t base_h; // base^h
    mpz_init( base_h );
    mpz_pow_ui( base_h, base, h );
    mpz_t base_raised; // base_h^i
    mpz_init( base_raised );
    for ( i = 0; i < D_count; i++ ) { // each iteration computes one term and then adds to B
      if ( i == 0 )
        mpz_set_ui( base_raised, 1 );
      else if ( i == 1 )
        mpz_set( base_raised, base_h );
        mpz_mul( base_raised, base_raised, base_h );
      mpz_mul( tempz1, base_raised, D[i] );
      mpz_add( B, B, tempz1 );
    mpz_clear( base_raised );
    mpz_clear( base_h );
  if ( mpz_cmp( base, sqrfreebase ) != 0 ) {
    mpz_divexact( tempz1, base, sqrfreebase );
    mpz_sqrt( tempz1, tempz1 );
    mpz_pow_ui( tempz1, tempz1, h );
    mpz_mul( B, B, tempz1 );
  mpz_pow_ui( tempz1, sqrfreebase, k );
  mpz_mul( B, B, tempz1 );
  mpz_sub( L, A, B );
  mpz_add( M, A, B );
  mpz_clear( B );
  mpz_clear( A );
  mpz_clear( tempz1 );
  for ( i = 0; i < C_count; i++ )
    mpz_clear( C[i] );
  free( C );
  C = NULL;
  for ( i = 0; i < D_count; i++ )
    mpz_clear( D[i] );
  free( D );
  D = NULL;
  return 0;
char quickprimecheck( mpz_t the_number ) {
  char retval = 'C';
  if ( mpz_cmp_ui( the_number, 1 ) == 0 ) {
    retval = 'N';
    return retval;
  int factor_type = mpz_probab_prime_p( the_number, 30 );
  if ( factor_type != 0 )
    retval = 'P';
  return retval;
// divide out denominator from running_N
void TFDivideOut( mpz_t running_N, long ldenominator, char* running_N_status, mpz_t square_root, struct factor_infos* Factor_Infos ) {
  mpz_t denominator;
  mpz_init_set_ui( denominator, ldenominator );
  long occurrences = ComputeOccurrences( running_N, denominator );
  long i = 0;
  for ( i = occurrences; i > 0; i-- )
    mpz_divexact_ui( running_N, running_N, ldenominator );
  *running_N_status = quickprimecheck( running_N );
  if ( *running_N_status == 'N' )
    mpz_set_ui( square_root, 1 );
    mpz_sqrt( square_root, running_N );
  AddFactorInfo( Factor_Infos, denominator, occurrences, 'P', NumberOfDigits( denominator ) );
  mpz_clear( denominator );
// Compute the prime factorization of a general number where square_root( the_number ) < MAXLONG - 255
void Factorize( mpz_t the_number, struct factor_infos* Factor_Infos ) {
  if ( Factor_Infos == NULL )
    return;
  // Wheel Factorization. eg. http://programmingpraxis.com/2009/05/08/wheel-factorization/
  mpz_t running_N;
  mpz_init_set( running_N, the_number );
  char running_N_status = quickprimecheck( running_N );
  if ( running_N_status == 'P' ) {
    AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
    mpz_clear( running_N );
    return;
  mpz_t square_root;
  mpz_init( square_root );
  mpz_sqrt( square_root, running_N );
  if ( mpz_cmp_ui( square_root, MAXLONG - 255 ) > 0 ) {
    printf( "Internal Error:  Trying to factorize too big a number.\n" );
    mpz_clear( square_root );
    mpz_clear( running_N );
    return;
  if ( mpz_divisible_ui_p( running_N, 2 ) != 0 ) {
    TFDivideOut( running_N, 2, &running_N_status, square_root, Factor_Infos );
    if ( running_N_status == 'P' ) {
      AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
      mpz_clear( square_root );
      mpz_clear( running_N );
      return;
  if ( mpz_divisible_ui_p( running_N, 3 ) != 0 ) {
    TFDivideOut( running_N, 3, &running_N_status, square_root, Factor_Infos );
    if ( running_N_status == 'P' ) {
      AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
      mpz_clear( square_root );
      mpz_clear( running_N );
      return;
  if ( mpz_divisible_ui_p( running_N, 5 ) != 0 ) {
    TFDivideOut( running_N, 5, &running_N_status, square_root, Factor_Infos );
    if ( running_N_status == 'P' ) {
      AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
      mpz_clear( square_root );
      mpz_clear( running_N );
      return;
  if ( mpz_divisible_ui_p( running_N, 7 ) != 0 ) {
    TFDivideOut( running_N, 7, &running_N_status, square_root, Factor_Infos );
    if ( running_N_status == 'P' ) {
      AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
      mpz_clear( square_root );
      mpz_clear( running_N );
      return;
  unsigned long i = 0;
  unsigned long denominator = 11;
        running_N_status == 'C';
        denominator += wheel7[i], i = ( i == 47 ? 0 : i + 1 ) ) {
      if ( mpz_divisible_ui_p( running_N, denominator ) != 0 )
        TFDivideOut( running_N, denominator, &running_N_status, square_root, Factor_Infos );
  if ( running_N_status == 'P' )
    AddFactorInfo( Factor_Infos, running_N, 1, 'P', NumberOfDigits( running_N ) );
  mpz_clear( square_root );
  mpz_clear( running_N );
void PermuteFactors( struct factor_infos* permutefactors, struct factor_infos* allfactors, long lefttumbler, mpz_t leftproduct, mpz_t the_number ) {
  long tumbler = lefttumbler + 1;
  mpz_t currfactor;
  mpz_init( currfactor );
  mpz_t currproduct;
  mpz_init( currproduct );
  long exponent;
  for ( exponent = 0; exponent <= permutefactors->the_factors[tumbler].occurrences; exponent++ ) {
    // compute primefactor^exponent
    if ( exponent == 0 )
      mpz_set_ui( currfactor, 1 );
    else if ( exponent == 1 )
      mpz_set( currfactor, permutefactors->the_factors[tumbler].the_factor );
      mpz_mul( currfactor, currfactor, permutefactors->the_factors[tumbler].the_factor );
    mpz_mul( currproduct, leftproduct, currfactor );
    if ( tumbler < permutefactors->count - 1 )
      PermuteFactors( permutefactors, allfactors, tumbler, currproduct, the_number );
    else if ( mpz_cmp_ui( currproduct, 1 ) != 0 && mpz_cmp( currproduct, the_number ) != 0 ) {
      char status = quickprimecheck( currproduct );
      long occurrences = ComputeOccurrences( the_number, currproduct );
      AddFactorInfo( allfactors, currproduct, occurrences, status, NumberOfDigits( currproduct ) );
  mpz_clear( currproduct );
  mpz_clear( currfactor );
// Compute all factors of the_number (Trivial factors of 1 and N defaulting to not be included)
void ComputeAllFactors( mpz_t N, struct factor_infos* Factor_Infos, int AddTrivialFactors ) {
  if ( Factor_Infos == NULL )
    return;
  if ( quickprimecheck( N ) == 'C' ) {
    struct factor_infos primes_only;
    Init_Factor_Infos( &primes_only );
    Factorize( N, &primes_only );
    mpz_t leftproduct;
    mpz_init_set_ui( leftproduct, 1 );
    PermuteFactors( &primes_only, Factor_Infos, -1, leftproduct, N );
    mpz_clear( leftproduct );
    Cleanup_Factor_Infos( &primes_only );
  if ( AddTrivialFactors ) {
    mpz_t tempZ;
    mpz_init_set_ui( tempZ, 1 );
    AddFactorInfo( Factor_Infos, tempZ, 1, 'N', 1 );
    if ( mpz_cmp_ui( N, 1 ) != 0 )
      AddFactorInfo( Factor_Infos, N, 1, '?', NumberOfDigits( N ) );
    mpz_clear( tempZ );
  qsort( Factor_Infos->the_factors, Factor_Infos->count, sizeof(struct factor_info), factor_info_cmpfunc );
// Compute the number of times denominator divides evenly into numerator
long ComputeOccurrences( mpz_t numerator, mpz_t denominator ) {
  // not handled values
  if ( mpz_cmp_ui( denominator, 1 ) <= 0 )
    return -1;
  long occurrences = 0;
  mpz_t quotient;
  mpz_init_set( quotient, numerator );
  while ( mpz_divisible_p( quotient, denominator ) ) {
    mpz_divexact( quotient, quotient, denominator );
    occurrences++;
  mpz_clear( quotient );
  return occurrences;
// Compute the Mobius function
int Mobius( mpz_t n ) {
  if ( mpz_cmp_ui( n, 1 ) == 0 )
    return 1;
  long total_occurrences = 0;
  int has_duplicate_occurrences = 0;
  if ( quickprimecheck( n ) == 'P' ) {
    total_occurrences++;
    struct factor_infos Factor_Infos;
    Init_Factor_Infos( &Factor_Infos );
View remainder of file in raw view
Provide feedback

Saved searches

Use saved searches to filter your results more quickly

FilesExpand file tree

cyclop.cpp

Latest commit

History

cyclop.cpp

File metadata and controls
