Added the Moebius function

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\_}

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);

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
 
 /*

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 */

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;

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

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
 };