Roots of unity: Difference between revisions

From Rosetta Code
Content added Content deleted
(added Fortran)
Line 256: Line 256:
=={{header|Python}}==
=={{header|Python}}==
{{works with|Python|2.5}}
{{works with|Python|2.5}}
<python>
import cmath


# format complex to 4 decimal places
Function <code>nthroots()</code> returns all n-th roots of unity.
class Complex(complex):
<python>from cmath import exp, pi
def nthroots(n):
def __repr__(self):
re_fmt = "%7.4f"
return [exp(k * 2j * pi / n) for k in range(1, n + 1)]</python>
im_fmt = ("+" if self.imag > 0 else "-") + "%6.4fj"
Example:
return (re_fmt+im_fmt) % (self.real, abs(self.imag))
>>> f = nthroots

>>> for n in range(1, 4):
for root in range(2, 10+1):
... print map(lambda c: ("%.3f " + ("+" if c.imag > 0 else "-") + " %.3gj") % (c.real, abs(c.imag)), f(n))
print "%3d"%root, [Complex(cmath.exp(2j*cmath.pi*n/root)) for n in range(1, root/2+1)]
['1.000 - 2.45e-016j']
</python>
['-1.000 + 1.22e-016j', '1.000 - 2.45e-016j']
Output:
['-0.500 + 0.866j', '-0.500 - 0.866j', '1.000 - 2.45e-016j']
2 [-1.0000+0.0000j]
3 [-0.5000+0.8660j]
4 [ 0.0000+1.0000j, -1.0000+0.0000j]
5 [ 0.3090+0.9511j, -0.8090+0.5878j]
6 [ 0.5000+0.8660j, -0.5000+0.8660j, -1.0000+0.0000j]
7 [ 0.6235+0.7818j, -0.2225+0.9749j, -0.9010+0.4339j]
8 [ 0.7071+0.7071j, 0.0000+1.0000j, -0.7071+0.7071j, -1.0000+0.0000j]
9 [ 0.7660+0.6428j, 0.1736+0.9848j, -0.5000+0.8660j, -0.9397+0.3420j]
10 [ 0.8090+0.5878j, 0.3090+0.9511j, -0.3090+0.9511j, -0.8090+0.5878j, -1.0000+0.0000j]


=={{header|Seed7}}==
=={{header|Seed7}}==

Revision as of 12:41, 9 August 2008

Task
Roots of unity
You are encouraged to solve this task according to the task description, using any language you may know.

The purpose of this task is to explore working with complex numbers. Given n, find the n-th roots of unity.

Ada

with Ada.Text_Io; use Ada.Text_Io;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Elementary_Functions;
with Ada.Text_Io.Complex_Io;

procedure Roots_Of_Unity is
   type Real is digits 10;
   package Complex_Type is new Ada.Numerics.Generic_Complex_Types(Real);
   use Complex_Type;
   package Complex_Functions is new Ada.Numerics.Generic_Complex_Elementary_Functions(Complex_Type);
   use Complex_Functions;
   package Complex_Io is new Ada.Text_Io.Complex_Io(Complex_Type);
   Factor : Complex;
begin
   for Root in 2..10 loop
      Put(Item => Root, Width => 4);
      Put(" ");
      for K in 1..Root/2 loop
         Factor := Exp(Real(K) * 2.0 * Ada.Numerics.Pi / Real(Root));
         Complex_Io.Put(Item => Factor, Fore => 2, Aft  => 4, Exp  => 0);
      end loop;
      New_Line;
   end loop;
end Roots_Of_Unity;

Output: 

  2 (-1.0000, 0.0000)
  3 (-0.5000, 0.8660)
  4 ( 0.0000, 1.0000)(-1.0000, 0.0000)
  5 ( 0.3090, 0.9511)(-0.8090, 0.5878)
  6 ( 0.5000, 0.8660)(-0.5000, 0.8660)(-1.0000, 0.0000)
  7 ( 0.6235, 0.7818)(-0.2225, 0.9749)(-0.9010, 0.4339)
  8 ( 0.7071, 0.7071)( 0.0000, 1.0000)(-0.7071, 0.7071)(-1.0000, 0.0000)
  9 ( 0.7660, 0.6428)( 0.1736, 0.9848)(-0.5000, 0.8660)(-0.9397, 0.3420)
 10 ( 0.8090, 0.5878)( 0.3090, 0.9511)(-0.3090, 0.9511)(-0.8090, 0.5878)(-1.0000, 0.0000)

ALGOL 68

FORMAT complex fmt=$g(-6,4)"⊥"g(-6,4)$;
FOR root FROM 2 TO 10 DO
  printf(($g(4)$,root));
  FOR n TO root OVER 2 DO
    printf(($xf(complex fmt)$,complex exp( 0 I 2*pi*n/root)))
  OD;
  printf($l$)
OD

Output: 

 +2 -1.000⊥0.0000
 +3 -.5000⊥0.8660
 +4 0.0000⊥1.0000 -1.000⊥0.0000
 +5 0.3090⊥0.9511 -.8090⊥0.5878
 +6 0.5000⊥0.8660 -.5000⊥0.8660 -1.000⊥0.0000
 +7 0.6235⊥0.7818 -.2225⊥0.9749 -.9010⊥0.4339
 +8 0.7071⊥0.7071 0.0000⊥1.0000 -.7071⊥0.7071 -1.000⊥0.0000
 +9 0.7660⊥0.6428 0.1736⊥0.9848 -.5000⊥0.8660 -.9397⊥0.3420
+10 0.8090⊥0.5878 0.3090⊥0.9511 -.3090⊥0.9511 -.8090⊥0.5878 -1.000⊥0.0000

BASIC

Works with: QuickBasic version 4.5
Translation of: Java

For high n's, this may repeat the root of 1 + 0*i. 

CLS
PI = 3.1415926#
n = 5 'this can be changed for any desired n
angle = 0 'start at angle 0
DO
	real = COS(angle) 'real axis is the x axis
	IF (ABS(real) < 10 ^ -5) THEN real = 0 'get rid of annoying sci notation
	imag = SIN(angle) 'imaginary axis is the y axis
	IF (ABS(imag) < 10 ^ -5) THEN imag = 0 'get rid of annoying sci notation
	PRINT real; "+"; imag; "i" 'answer on every line
	angle = angle + (2 * PI) / n
'all the way around the circle at even intervals
LOOP WHILE angle < 2 * PI

C++

<cpp>#include <complex>

  1. include <cmath>
  2. include <iostream>

double const pi = 4 * std::atan(1);

int main() { 

 for (int n = 2; n <= 10; ++n)
 {
   std::cout << n << ": ";
   for (int k = 0; k < n; ++k)
     std::cout << std::polar(1, 2*pi*k/n) << " ";
   std::cout << std::endl;
 }

}</cpp>

D

Works with: D version 2.012
Works with: D version 1.028

<d>module nthroots ; import std.stdio, std.math ;

creal[] nthroots(int n) { 

 creal[] res ;
 for(int k = 1 ; k <= n ; k++)
   res ~= expi(PI*2*k/n) ;
 return res ;

} void main() { 

 for(int i = 1; i <= 8 ; i++)
   writefln("%2dth : %5.2f", i, nthroots(i)) ;

}</d>

Forth

Complex numbers are not a native type in Forth, so we calculate the roots by hand. 

: f0. ( f -- )
  fdup 0e 0.001e f~ if fdrop 0e then f. ;
: .roots ( n -- )
  dup 1 do
    pi i 2* 0 d>f f* dup 0 d>f f/          ( F: radians )
    fsincos cr ." real " f0. ." imag " f0.
  loop drop ;

3 set-precision
5 .roots

Fortran

Works with: Fortran version 90 and later
PROGRAM Roots

  COMPLEX :: root 
  INTEGER :: n
  REAL :: angle, pi

  pi = 4.0 * ATAN(1.0)
  DO n = 2, 7
    angle = 0.0
    WRITE(*,"(I1,A)", ADVANCE="NO") n,": "
    DO WHILE (angle < 6.283)
      root = CMPLX(COS(angle), SIN(angle))
      WRITE(*,"(SP,2F7.4,A)", ADVANCE="NO") root, "j  "
      angle = angle + (2.0*pi) / REAL(n)
    END DO
    WRITE(*,*)
  END DO

END PROGRAM Roots

Output 

2: +1.0000+0.0000j  -1.0000+0.0000j   
3: +1.0000+0.0000j  -0.5000+0.8660j  -0.5000-0.8660j   
4: +1.0000+0.0000j  +0.0000+1.0000j  -1.0000+0.0000j  +0.0000-1.0000j   
5: +1.0000+0.0000j  +0.3090+0.9511j  -0.8090+0.5878j  -0.8090-0.5878j  +0.3090-0.9511j   
6: +1.0000+0.0000j  +0.5000+0.8660j  -0.5000+0.8660j  -1.0000+0.0000j  -0.5000-0.8660j  +0.5000-0.8660j 
7: +1.0000+0.0000j  +0.6235+0.7818j  -0.2225+0.9749j  -0.9010+0.4339j  -0.9010-0.4339j  -0.2225-0.9749j  +0.6235-0.7818j

Haskell

import Data.Complex

rootsOfUnity n = [mkPolar 1.0 (2*pi*k/n) | k <- [1..n]]

Output: 

*Main> rootsOfUnity 3
[(-0.4999999999999998) :+ 0.8660254037844387,
 (-0.5000000000000004) :+ (-0.8660254037844384),
 1.0 :+ (-2.4492127076447545e-16)]

IDL

For some example n: 

 n = 5
 print,  exp( dcomplex( 0, 2*!pi/n) ) ^ ( 1 + indgen(n) )

Outputs: 

 ( 0.30901696, 0.95105653)( -0.80901704, 0.58778520)( -0.80901693, -0.58778534)( 0.30901713, -0.95105647)( 1.0000000, 1.7484556e-007)

J

   rou=: [: ^ i. * (o. 0j2) % ]

   rou 4
1 0j1 _1 0j_1

   rou 5
1 0.309017j0.951057 _0.809017j0.587785 _0.809017j_0.587785 0.309017j_0.951057

The computation can also be written as a loop, shown here for comparison only. 

rou1=: 3 : 0
 z=. 0 $ r=. ^ o. 0j2 % y [ e=. 1
 for. i.y do.
  z=. z,e
  e=. e*r
 end.
 z
)

Java

Java doesn't have a nice way of dealing with complex numbers, so the real and imaginary parts are calculated separately based on the angle and printed together. There are also checks in this implementation to get rid of extremely small values (< 1.0E-3 where scientific notation sets in for Doubles). Instead, they are simply represented as 0. To remove those checks (for very high n's), remove both if statements. <java>public static void unity(int n){ //all the way around the circle at even intervals for(double angle = 0;angle < 2 * Math.PI;angle += (2 * Math.PI) / n){ double real = Math.cos(angle); //real axis is the x axis if(Math.abs(real) < 1.0E-3) real = 0.0; //get rid of annoying sci notation double imag = Math.sin(angle); //imaginary axis is the y axis if(Math.abs(imag) < 1.0E-3) imag = 0.0; //get rid of annoying sci notation System.out.print(real + " + " + imag + "i\t"); //tab-separated answers } }</java>

OCaml

<ocaml>open Complex

let pi = 4. *. atan 1.

let () = 

 for n = 1 to 10 do
   Printf.printf "%2d " n;
   for k = 1 to n do
     let ret = polar 1. (2. *. pi *. float_of_int k /. float_of_int n) in
       Printf.printf "(%f + %f i)" ret.re ret.im
   done;
   print_newline ()
 done</ocaml>

Perl

Works with: Perl version 5.8.8
Library: Math::ComplexComplex

<perl>use Math::Complex;

foreach $n (1 .. 10) { 

 printf "%2d ", $n;
 foreach $k (1 .. $n) {
   $ret = cplxe(1, 2 * pi * $k / $n);
   $ret->display_format('style' => 'cartesian', 'format' => '%.3f');
   print "($ret)";
 }
 print "\n";

}</perl>

Output: 

 1 (1.000-0.000i)
 2 (-1.000+0.000i)(1.000-0.000i)
 3 (-0.500+0.866i)(-0.500-0.866i)(1.000-0.000i)
 4 (0.000+1.000i)(-1.000+0.000i)(-0.000-1.000i)(1.000-0.000i)
 5 (0.309+0.951i)(-0.809+0.588i)(-0.809-0.588i)(0.309-0.951i)(1.000-0.000i)
 6 (0.500+0.866i)(-0.500+0.866i)(-1.000+0.000i)(-0.500-0.866i)(0.500-0.866i)(1.000-0.000i)
 7 (0.623+0.782i)(-0.223+0.975i)(-0.901+0.434i)(-0.901-0.434i)(-0.223-0.975i)(0.623-0.782i)(1.000-0.000i)
 8 (0.707+0.707i)(0.000+1.000i)(-0.707+0.707i)(-1.000+0.000i)(-0.707-0.707i)(-0.000-1.000i)(0.707-0.707i)(1.000-0.000i)
 9 (0.766+0.643i)(0.174+0.985i)(-0.500+0.866i)(-0.940+0.342i)(-0.940-0.342i)(-0.500-0.866i)(0.174-0.985i)(0.766-0.643i)(1.000-0.000i)
10 (0.809+0.588i)(0.309+0.951i)(-0.309+0.951i)(-0.809+0.588i)(-1.000+0.000i)(-0.809-0.588i)(-0.309-0.951i)(0.309-0.951i)(0.809-0.588i)(1.000-0.000i)

Python

Works with: Python version 2.5

<python> import cmath

  1. format complex to 4 decimal places

class Complex(complex): 

 def __repr__(self):
   re_fmt = "%7.4f"
   im_fmt = ("+" if self.imag > 0 else "-") + "%6.4fj"
   return (re_fmt+im_fmt) % (self.real, abs(self.imag))

for root in range(2, 10+1): 

 print "%3d"%root, [Complex(cmath.exp(2j*cmath.pi*n/root)) for n in range(1, root/2+1)]

</python> Output: 

 2 [-1.0000+0.0000j]
 3 [-0.5000+0.8660j]
 4 [ 0.0000+1.0000j, -1.0000+0.0000j]
 5 [ 0.3090+0.9511j, -0.8090+0.5878j]
 6 [ 0.5000+0.8660j, -0.5000+0.8660j, -1.0000+0.0000j]
 7 [ 0.6235+0.7818j, -0.2225+0.9749j, -0.9010+0.4339j]
 8 [ 0.7071+0.7071j,  0.0000+1.0000j, -0.7071+0.7071j, -1.0000+0.0000j]
 9 [ 0.7660+0.6428j,  0.1736+0.9848j, -0.5000+0.8660j, -0.9397+0.3420j]
10 [ 0.8090+0.5878j,  0.3090+0.9511j, -0.3090+0.9511j, -0.8090+0.5878j, -1.0000+0.0000j]

Seed7

$ include "seed7_05.s7i";
  include "float.s7i";
  include "complex.s7i";

const proc: main is func
  local
    var integer: n is 0;
    var integer: k is 0;
  begin
    for n range 2 to 10 do
      write(n lpad 2 <& ": ");
      for k range 0 to pred(n) do
        write(polar(1.0, 2.0 * PI * flt(k) / flt(n)) digits 4 lpad 15 <& " ");
      end for;
      writeln;
    end for;
  end func;

Output: 

 2:  1.0000+0.0000i -1.0000+0.0000i
 3:  1.0000+0.0000i -0.5000+0.8660i -0.5000-0.8660i
 4:  1.0000+0.0000i  0.0000+1.0000i -1.0000+0.0000i  0.0000-1.0000i
 5:  1.0000+0.0000i  0.3090+0.9511i -0.8090+0.5878i -0.8090-0.5878i  0.3090-0.9511i
 6:  1.0000+0.0000i  0.5000+0.8660i -0.5000+0.8660i -1.0000+0.0000i -0.5000-0.8660i  0.5000-0.8660i
 7:  1.0000+0.0000i  0.6235+0.7818i -0.2225+0.9749i -0.9010+0.4339i -0.9010-0.4339i -0.2225-0.9749i  0.6235-0.7818i
 8:  1.0000+0.0000i  0.7071+0.7071i  0.0000+1.0000i -0.7071+0.7071i -1.0000+0.0000i -0.7071-0.7071i  0.0000-1.0000i  0.7071-0.7071i
 9:  1.0000+0.0000i  0.7660+0.6428i  0.1736+0.9848i -0.5000+0.8660i -0.9397+0.3420i -0.9397-0.3420i -0.5000-0.8660i  0.1736-0.9848i  0.7660-0.6428i
10:  1.0000+0.0000i  0.8090+0.5878i  0.3090+0.9511i -0.3090+0.9511i -0.8090+0.5878i -1.0000+0.0000i -0.8090-0.5878i -0.3090-0.9511i  0.3090-0.9511i  0.8090-0.5878i
Retrieved from "https://rosettacode.org/wiki/Roots_of_unity?oldid=50857"