diff --git a/doc/manual/functions.tex b/doc/manual/functions.tex index 566c2fab..606ce491 100644 --- a/doc/manual/functions.tex +++ b/doc/manual/functions.tex @@ -735,6 +735,30 @@ \section{mod2\_}\index{mod2\_}\index{function!mod2\_} $+[(p+1)/2]$. %--#] mod2_ : +%--#[ moebius_ : + +\section{moebius\_}\index{moebius\_}\index{function!moebius\_} +\label{funmoebius} +\noindent The Moebius function. It has the property: +\begin{itemize} +\item moebius\_(1) = 1 +\item moebius\_(n) = 0 if n can be divided by the square of a prime number +\item moebius\_(n) = 1 if n contains an even number of different prime factors +\item moebius\_(n) = -1 if n contains an odd number of different prime factors +\end{itemize} +The number n must be a positive integer for the evaluation to take place. +The argument should be no bigger than can be allowed by the biggest primes +in Form. This means for a 64-bit computer arguments above $2^{31}$ may not +give a correct result with the current algorithms that try for prime +factors. There is a second algorithm that does not suffer from this +problem: +\begin{verbatim} + id moebius_(n) = 1-sum_(j,1,n-1,integer_(n/j)*moebius_(j)); +\end{verbatim} +but this algorithm needs all lower values and hence is inherently +quadratic. Therefore we did not use it. + +%--#] moebius_ : %--#[ mul_ : \section{mul\_}\index{mul\_}\index{function!mul\_} diff --git a/sources/declare.h b/sources/declare.h index 0cbf47be..175c9e60 100644 --- a/sources/declare.h +++ b/sources/declare.h @@ -1475,6 +1475,7 @@ extern int TestFunFlag(PHEAD WORD *); extern WORD CompareSymbols(WORD *,WORD *,WORD); extern WORD CompareHSymbols(WORD *,WORD *,WORD); extern WORD NextPrime(PHEAD WORD); +extern WORD Moebius(PHEAD WORD); extern UWORD wranf(PHEAD0); extern UWORD iranf(PHEAD UWORD); extern void iniwranf(PHEAD0); diff --git a/sources/ftypes.h b/sources/ftypes.h index d2dcb427..7637dfae 100644 --- a/sources/ftypes.h +++ b/sources/ftypes.h @@ -449,14 +449,15 @@ typedef int (*TFUN1)(); #define PERMUTATIONS 103 #define PARTITIONS 104 #define MULFUNCTION 105 -#define TOPOLOGIES 106 -#define DIAGRAMS 107 -#define VERTEX 108 -#define EDGE 109 +#define MOEBIUS 106 +#define TOPOLOGIES 107 +#define DIAGRAMS 108 +#define VERTEX 109 +#define EDGE 110 /*#define ALLWILD 109 ???? */ -#define SIZEOFFUNCTION 110 +#define SIZEOFFUNCTION 111 -#define MAXBUILTINFUNCTION 110 +#define MAXBUILTINFUNCTION 111 #define FIRSTUSERFUNCTION 150 /* diff --git a/sources/inivar.h b/sources/inivar.h index e2676507..7c7d3bf0 100644 --- a/sources/inivar.h +++ b/sources/inivar.h @@ -231,6 +231,7 @@ static struct fixedfun { ,{"perm_" ,1 ,0 ,0 ,0} /* PERMUTATIONS */ ,{"partitions_" ,1 ,0 ,0 ,0} /* PARTITIONS */ ,{"mul_" ,0 ,0 ,0 ,0} /* MULFUNCTION */ + ,{"moebius_" ,0 ,0 ,0 ,0} /* MOEBIUS */ ,{"topologies_" ,0 ,0 ,0 ,0} /* TOPOLOGIES */ ,{"diagrams_" ,0 ,0 ,0 ,0} /* DIAGRAMS */ ,{"node_" ,0 ,0 ,0 ,0} /* VERTEX */ diff --git a/sources/normal.c b/sources/normal.c index 0b875f67..b6fb9b9e 100644 --- a/sources/normal.c +++ b/sources/normal.c @@ -1959,6 +1959,20 @@ nospec: pcom[ncom++] = t; } else goto defaultcase; break; + case MOEBIUS: +/* + Numerical function giving -1,0,1 or no evaluation +*/ + if ( t[1] == FUNHEAD+2 && t[FUNHEAD] == -SNUMBER + && t[FUNHEAD+1] > 0 ) { + WORD val = Moebius(BHEAD t[FUNHEAD+1]); + if ( val == 0 ) goto NormZero; + if ( val < 0 ) ncoef = -ncoef; + } + else { + pcom[ncom++] = t; + } + break; case LNUMBER : ncoef = REDLENG(ncoef); if ( Mully(BHEAD (UWORD *)n_coef,&ncoef,(UWORD *)(t+3),t[2]) ) goto FromNorm; diff --git a/sources/reken.c b/sources/reken.c index 8b72ed5d..88406883 100644 --- a/sources/reken.c +++ b/sources/reken.c @@ -3708,6 +3708,66 @@ nexti:; /* #] NextPrime : + #[ Moebius : + + Returns the value of the Moebius function if n fits inside + a WORD. + The method we use is a bit like the sieve of Erathostenes. +*/ + +WORD Moebius(PHEAD WORD nn) +{ + WORD i,n = nn, x; + LONG newsize; + char *newtable, mu; +#if ( BITSINWORD == 32 ) + if ( AR.notfirstprime == 0 ) StartPrimeList(BHEAD0); +#endif +/* + First we make sure that: + a: the table is big enough. + b: the number is not already in the table. +*/ + if ( nn >= AR.moebiustablesize ) { + if ( AR.moebiustablesize <= 0 ) { newsize = nn + 20; } + else { newsize = nn*2; } + if ( newsize > MAXPOSITIVE ) newsize = MAXPOSITIVE; + newtable = (char *)Malloc1(newsize*sizeof(char),"Moebius"); + for ( i = 0; i < AR.moebiustablesize; i++ ) newtable[i] = AR.moebiustable[i]; + for ( ; i < newsize; i++ ) newtable[i] = 2; + if ( AR.moebiustablesize > 0 ) M_free(AR.moebiustable,"Moebius"); + AR.moebiustable = newtable; + AR.moebiustablesize = newsize; + } + if ( AR.moebiustable[nn] != 2 ) return((WORD)AR.moebiustable[nn]); + mu = 1; + if ( n == 1 ) goto putvalue; + if ( n % 2 == 0 ) { + n /= 2; + if ( n % 2 == 0 ) { mu = 0; goto putvalue; } + if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n]; goto putvalue; } + mu = -mu; + if ( n == 1 ) goto putvalue; + } + for ( i = 0; i < AR.numinprimelist; i++ ) { + x = AR.PrimeList[i]; + if ( n % x == 0 ) { + n /= x; + if ( n % x == 0 ) { mu = 0; goto putvalue; } + if ( AR.moebiustable[n] != 2 ) { mu = -AR.moebiustable[n]; goto putvalue; } + mu = -mu; + if ( n == 1 ) goto putvalue; + } + if ( n < x*x ) break; /* notice that x*x always fits inside a WORD */ + } + mu = -mu; +putvalue: + AR.moebiustable[nn] = mu; + return((WORD)mu); +} + +/* + #] Moebius : #[ wranf : A random number generator that generates random WORDs with a very diff --git a/sources/structs.h b/sources/structs.h index 44c9df1e..401bcd97 100644 --- a/sources/structs.h +++ b/sources/structs.h @@ -1926,6 +1926,7 @@ struct R_const { WORD *CompressPointer; /* (R) */ COMPARE CompareRoutine; ULONG *wranfia; + char *moebiustable; LONG OldTime; /* (R) Zero time. Needed in timer. */ LONG InInBuf; /* (R) Characters in input buffer. Scratch files. */ @@ -1973,17 +1974,18 @@ struct R_const { WORD CurExpr; /* (S) Number of current expression */ WORD SortType; /* A copy of AC.SortType to play with */ WORD ShortSortCount; /* For On FewerStatistics 10; */ + WORD moebiustablesize; #if ( BITSINWORD == 32 ) #ifdef WITHPTHREADS - PADPOSITION(8,7,8,5026,0); + PADPOSITION(8,7,8,5027,0); #else - PADPOSITION(8,7,7,5026,0); + PADPOSITION(8,7,7,5027,0); #endif #else #ifdef WITHPTHREADS - PADPOSITION(8,7,8,24,0); + PADPOSITION(8,7,8,25,0); #else - PADPOSITION(8,7,7,24,0); + PADPOSITION(8,7,7,25,0); #endif #endif };