Möbius function: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎{{header|REXX}}: added a REXX stub.)
(→‎{{header|REXX}}: added the REXX computer programming language for this task.)
Line 161: Line 161:


=={{header|REXX}}==
=={{header|REXX}}==
<lang rexx>/*REXX pgm computes and shows a value grid of the Möbius funtion for a range of integers*/
<lang rexx></lang>
parse arg LO HI grp . /*obtain optional arguments from the CL*/
if LO=='' | LO=="," then LO= 0 /*Not specified? Then use the default.*/
if HI=='' | HI=="," then HI= 199 /* " " " " " " */
if grp=='' | grp=="," then grp= 20 /* " " " " " " */
/* ______ */
call genP HI /*generate primes up to the √ HI */
say center(' The Möbius sequence from ' LO " ──► " HI" ", max(50, grp*3), '═') /*title*/
$= /*variable holds output grid of GRP #s.*/
do j=LO to HI; $= $ right(mobius(j), 2) /*process some numbers from LO ──► HI.*/
if words($)==grp then do; say substr($, 2); $= /*show grid if fully populated,*/
end /* and nullify it for more #s.*/
end /*j*/ /*for small grids, using wordCnt is OK.*/

if $\=='' then say substr($, 2) /*handle any residual numbers not shown*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
mobius: procedure expose @.; parse arg x /*obtain a integer to be tested for mu.*/
if x==0 then return '∙' /*Zero? Then return a symbol for null.*/
if x==1 then return 1 /*Unity? Then return unity. */
#= 0 /*start with a value of zero. */
do k=1; p= @.k /*get the Kth (pre─generated) prime.*/
if p>x then leave /*prime (P) > X? Then we're done. */
if p*p>x then do; #= #+1; leave /*prime (P**2 > X? Bump # and leave.*/
end
if x//p==0 then do; #= #+1 /*X divisible by P? Bump mu number. */
x= x % p /* Divide by prime. */
if x//p==0 then return 0 /*X÷by P? Then return zero.*/
end
end /*k*/
return -1 ** # /*raise negative unity to the mu power.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP: @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6= 13; nP=6 /*assign low primes; # primes. */
do lim=nP until lim*lim>=HI; end /*only keep primes up to the sqrt(HI). */
do j=@.nP+4 by 2 to lim /*only find odd primes from here on. */
if j// 3==0 then iterate /*is J divisible by #3 Then not prime*/
parse var j '' -1 _;if _==5 then iterate /*Is last digit a "5"? " " " */
if j// 7==0 then iterate /*is J divisible by 7? " " " */
if j//11==0 then iterate /* " " " " 11? " " " */
if j//13==0 then iterate /*is " " " 13? " " " */
do k=7 while k*k<=j /*divide by some generated odd primes. */
if j // @.k==0 then iterate j /*Is J divisible by P? Then not prime*/
end /*k*/ /* [↓] a prime (J) has been found. */
nP= nP+1; if nP<=lim then @.nP=j /*bump prime count; assign prime to @.*/
end /*j*/; return</lang>
{{out|output|text=&nbsp; when using the default inputs:}}

Output note: &nbsp; note the use of a bullet (&bull;) to signify that a "null" is being shown (for the 0<sup>th</sup> entry).
<pre>
══════════ The Möbius sequence from 0 ──► 199 ═══════════
∙ 1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1
0 1 1 -1 0 0 1 0 0 -1 -1 -1 0 1 1 1 0 -1 1 1
0 -1 -1 -1 0 0 1 -1 0 0 0 1 0 -1 0 1 0 1 1 -1
0 -1 1 0 0 1 -1 -1 0 1 -1 -1 0 -1 1 0 0 1 -1 -1
0 0 1 -1 0 1 1 1 0 -1 0 1 0 1 1 1 0 -1 0 0
0 -1 -1 -1 0 -1 1 -1 0 -1 -1 1 0 -1 -1 1 0 0 1 1
0 0 1 1 0 0 0 -1 0 1 -1 -1 0 1 1 0 0 -1 -1 -1
0 1 1 1 0 1 1 0 0 -1 0 -1 0 0 -1 1 0 -1 1 1
0 1 0 -1 0 -1 1 -1 0 0 -1 0 0 1 -1 0 0 1 1 1
0 1 -1 1 0 1 -1 1 0 0 -1 1 0 1 1 -1 0 1 0 1
</pre>

Revision as of 21:34, 25 January 2020

Möbius function is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The classical Möbius function: μ(n) is an important multiplicative function in number theory and combinatorics.

There are several ways to implement a Möbius function.

A fairly straightforward method is to find the prime factors of a positive integer n, then define μ(n) based on the sum of the primitive factors. It has the values {−1, 0, 1} depending on the factorization of n:

  • μ(1) is defined to be 1.
  • μ(n) = 1 if n is a square-free positive integer with an even number of prime factors.
  • μ(n) = −1 if n is a square-free positive integer with an odd number of prime factors.
  • μ(n) = 0 if n has a squared prime factor.


Task
  • Write a routine (function, procedure, whatever) μ(n) to find the Möbius number for a positive integer n.
  • Use that routine to find and display here, on this page, at least the first 99 terms in a grid layout. (Not just one long line or column of numbers.)


See also


Related Tasks


Factor

The mobius word exists in the math.extras vocabulary. See the implementation here.

Works with: Factor version 0.99 2020-01-23

<lang factor>USING: formatting grouping io math.extras math.ranges sequences ;

"First 199 terms of the Möbius sequence:" print 199 [1,b] [ mobius ] map " " prefix 20 group [ [ "%3s" printf ] each nl ] each</lang>

Output:
First 199 terms of the Möbius sequence:
     1 -1 -1  0 -1  1 -1  0  0  1 -1  0 -1  1  1  0 -1  0 -1
  0  1  1 -1  0  0  1  0  0 -1 -1 -1  0  1  1  1  0 -1  1  1
  0 -1 -1 -1  0  0  1 -1  0  0  0  1  0 -1  0  1  0  1  1 -1
  0 -1  1  0  0  1 -1 -1  0  1 -1 -1  0 -1  1  0  0  1 -1 -1
  0  0  1 -1  0  1  1  1  0 -1  0  1  0  1  1  1  0 -1  0  0
  0 -1 -1 -1  0 -1  1 -1  0 -1 -1  1  0 -1 -1  1  0  0  1  1
  0  0  1  1  0  0  0 -1  0  1 -1 -1  0  1  1  0  0 -1 -1 -1
  0  1  1  1  0  1  1  0  0 -1  0 -1  0  0 -1  1  0 -1  1  1
  0  1  0 -1  0 -1  1 -1  0  0 -1  0  0 -1 -1  0  0  1  1 -1
  0 -1 -1  1  0  1 -1  1  0  0 -1 -1  0 -1  1 -1  0 -1  0 -1

Go

<lang go>package main

import "fmt"

func möbius(to int) []int {

   if to < 1 {
       to = 1
   }
   mobs := make([]int, to+1) // all zero by default
   primes := []int{2}
   for i := 1; i <= to; i++ {
       j := i
       cp := 0      // counts prime factors
       spf := false // true if there is a square prime factor
       for _, p := range primes {
           if p > j {
               break
           }
           if j%p == 0 {
               j /= p
               cp++
           }
           if j%p == 0 {
               spf = true
               break
           }
       }
       if cp == 0 && i > 2 {
           cp = 1
           primes = append(primes, i)
       }
       if !spf {
           if cp%2 == 0 {
               mobs[i] = 1
           } else {
               mobs[i] = -1
           }
       }
   }
   return mobs

}

func main() {

   mobs := möbius(199)
   fmt.Println("Möbius sequence - First 199 terms:")
   for i := 0; i < 200; i++ {
       if i == 0 {
           fmt.Print("    ")
           continue
       }
       if i%20 == 0 {
           fmt.Println()
       }
       fmt.Printf("  % d", mobs[i])
   }

}</lang>

Output:
Möbius sequence - First 199 terms:
       1  -1  -1   0  -1   1  -1   0   0   1  -1   0  -1   1   1   0  -1   0  -1
   0   1   1  -1   0   0   1   0   0  -1  -1  -1   0   1   1   1   0  -1   1   1
   0  -1  -1  -1   0   0   1  -1   0   0   0   1   0  -1   0   1   0   1   1  -1
   0  -1   1   0   0   1  -1  -1   0   1  -1  -1   0  -1   1   0   0   1  -1  -1
   0   0   1  -1   0   1   1   1   0  -1   0   1   0   1   1   1   0  -1   0   0
   0  -1  -1  -1   0  -1   1  -1   0  -1  -1   1   0  -1  -1   1   0   0   1   1
   0   0   1   1   0   0   0  -1   0   1  -1  -1   0   1   1   0   0  -1  -1  -1
   0   1   1   1   0   1   1   0   0  -1   0  -1   0   0  -1   1   0  -1   1   1
   0   1   0  -1   0  -1   1  -1   0   0  -1   0   0  -1  -1   0   0   1   1  -1
   0  -1  -1   1   0   1  -1   1   0   0  -1  -1   0  -1   1  -1   0  -1   0  -1

Perl 6

Works with: Rakudo version 2019.11

Möbius number is not defined for n == 0. Perl 6 arrays are indexed from 0 so store a blank value at position zero to keep n and μ(n) aligned.

<lang perl6>use Prime::Factor;

sub μ (Int \n) {

   my @p = prime-factors(n);
   +@p == +Bag(@p).keys ?? +@p %% 2 ?? 1 !! -1 !! 0

}

my @möbius = lazy flat ' ', 1, (2..*).map: -> \n { μ(n) };

  1. The Task

put "Möbius sequence - First 199 terms:\n",

   @möbius[^200]».fmt('%3s').rotor(20).join("\n");</lang>
Output:
Möbius sequence - First 199 terms:
      1  -1  -1   0  -1   1  -1   0   0   1  -1   0  -1   1   1   0  -1   0  -1
  0   1   1  -1   0   0   1   0   0  -1  -1  -1   0   1   1   1   0  -1   1   1
  0  -1  -1  -1   0   0   1  -1   0   0   0   1   0  -1   0   1   0   1   1  -1
  0  -1   1   0   0   1  -1  -1   0   1  -1  -1   0  -1   1   0   0   1  -1  -1
  0   0   1  -1   0   1   1   1   0  -1   0   1   0   1   1   1   0  -1   0   0
  0  -1  -1  -1   0  -1   1  -1   0  -1  -1   1   0  -1  -1   1   0   0   1   1
  0   0   1   1   0   0   0  -1   0   1  -1  -1   0   1   1   0   0  -1  -1  -1
  0   1   1   1   0   1   1   0   0  -1   0  -1   0   0  -1   1   0  -1   1   1
  0   1   0  -1   0  -1   1  -1   0   0  -1   0   0  -1  -1   0   0   1   1  -1
  0  -1  -1   1   0   1  -1   1   0   0  -1  -1   0  -1   1  -1   0  -1   0  -1

REXX

<lang rexx>/*REXX pgm computes and shows a value grid of the Möbius funtion for a range of integers*/ parse arg LO HI grp . /*obtain optional arguments from the CL*/ if LO== | LO=="," then LO= 0 /*Not specified? Then use the default.*/ if HI== | HI=="," then HI= 199 /* " " " " " " */ if grp== | grp=="," then grp= 20 /* " " " " " " */

                                                /*                            ______   */

call genP HI /*generate primes up to the √ HI */ say center(' The Möbius sequence from ' LO " ──► " HI" ", max(50, grp*3), '═') /*title*/ $= /*variable holds output grid of GRP #s.*/

   do j=LO  to  HI;   $= $ right(mobius(j), 2)  /*process some numbers from  LO ──► HI.*/
   if words($)==grp  then do;  say substr($, 2);  $=    /*show grid if fully populated,*/
                          end                           /*  and nullify it for more #s.*/
   end   /*j*/                                  /*for small grids, using wordCnt is OK.*/

if $\== then say substr($, 2) /*handle any residual numbers not shown*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ mobius: procedure expose @.; parse arg x /*obtain a integer to be tested for mu.*/

       if x==0  then return '∙'                 /*Zero?  Then return a symbol for null.*/
       if x==1  then return  1                  /*Unity? Then return unity.            */
       #= 0                                     /*start with a value of zero.          */
          do k=1;  p= @.k                       /*get the  Kth  (pre─generated)  prime.*/
          if p>x  then leave                    /*prime (P)   > X?    Then we're done. */
          if p*p>x  then do;   #= #+1;  leave   /*prime (P**2 > X?    Bump # and leave.*/
                         end
          if x//p==0  then do; #= #+1           /*X divisible by P?   Bump mu number.  */
                               x= x % p         /*                    Divide by prime. */
                               if x//p==0  then return 0  /*X÷by P?   Then return zero.*/
                           end
          end   /*k*/
       return -1 ** #                           /*raise negative unity to the mu power.*/

/*──────────────────────────────────────────────────────────────────────────────────────*/ genP: @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6= 13; nP=6 /*assign low primes; # primes. */

      do lim=nP  until lim*lim>=HI;  end        /*only keep primes up to the sqrt(HI). */
      do j=@.nP+4  by 2  to lim                 /*only find odd primes from here on.   */
      if j// 3==0  then iterate                 /*is J divisible by  #3  Then not prime*/
      parse var j  -1 _;if _==5  then iterate /*Is last digit a "5"?     "   "    "  */
      if j// 7==0  then iterate                 /*is J divisible by  7?    "   "    "  */
      if j//11==0  then iterate                 /* " "     "      " 11?    "   "    "  */
      if j//13==0  then iterate                 /*is "     "      " 13?    "   "    "  */
         do k=7  while k*k<=j                   /*divide by some generated odd primes. */
         if j // @.k==0  then iterate j         /*Is J divisible by  P?  Then not prime*/
         end   /*k*/                            /* [↓]  a prime  (J)  has been found.  */
      nP= nP+1;  if nP<=lim  then @.nP=j        /*bump prime count; assign prime to  @.*/
      end      /*j*/;             return</lang>
output   when using the default inputs:

Output note:   note the use of a bullet (•) to signify that a "null" is being shown (for the 0th entry).

══════════ The Möbius sequence from  0  ──►  199 ═══════════
 ∙  1 -1 -1  0 -1  1 -1  0  0  1 -1  0 -1  1  1  0 -1  0 -1
 0  1  1 -1  0  0  1  0  0 -1 -1 -1  0  1  1  1  0 -1  1  1
 0 -1 -1 -1  0  0  1 -1  0  0  0  1  0 -1  0  1  0  1  1 -1
 0 -1  1  0  0  1 -1 -1  0  1 -1 -1  0 -1  1  0  0  1 -1 -1
 0  0  1 -1  0  1  1  1  0 -1  0  1  0  1  1  1  0 -1  0  0
 0 -1 -1 -1  0 -1  1 -1  0 -1 -1  1  0 -1 -1  1  0  0  1  1
 0  0  1  1  0  0  0 -1  0  1 -1 -1  0  1  1  0  0 -1 -1 -1
 0  1  1  1  0  1  1  0  0 -1  0 -1  0  0 -1  1  0 -1  1  1
 0  1  0 -1  0 -1  1 -1  0  0 -1  0  0  1 -1  0  0  1  1  1
 0  1 -1  1  0  1 -1  1  0  0 -1  1  0  1  1 -1  0  1  0  1