Sum of a series

From Rosetta Code
Jump to: navigation, search
Task
Sum of a series
You are encouraged to solve this task according to the task description, using any language you may know.
Compute the n-th term of a series, i.e. the sum of the n first terms of the corresponding sequence. Informally, it is also called the sum of the series, thus the title of this task.

For this task, use:

S_n = \sum_{k=1}^n \frac{1}{k^2}

and compute S1000.

This approximates the zeta function for s=2, whose exact value

\zeta(2) = {\pi^2\over 6}

is the solution of the Basel problem.

Contents

[edit] ACL2

(defun sum-x^-2 (max-x)
   (if (zp max-x)
       0
       (+ (/ (* max-x max-x))
          (sum-x^-2 (1- max-x)))))

[edit] ActionScript

function partialSum(n:uint):Number
{
	var sum:Number = 0;
	for(var i:uint = 1; i <= n; i++)
		sum += 1/(i*i);
	return sum;
}
trace(partialSum(1000));

[edit] Ada

with Ada.Text_Io; use Ada.Text_Io;
 
procedure Sum_Series is
   function F(X : Long_Float) return Long_Float is
   begin
      return 1.0 / X**2;
   end F;
   package Lf_Io is new Ada.Text_Io.Float_Io(Long_Float);
   use Lf_Io;
   Sum : Long_Float := 0.0;
   subtype Param_Range is Integer range 1..1000;
begin
   for I in Param_Range loop
      Sum := Sum + F(Long_Float(I));
   end loop;
   Put("Sum of F(x) from" & Integer'Image(Param_Range'First) &
      " to" & Integer'Image(Param_Range'Last) & " is ");
   Put(Item => Sum, Aft => 10, Exp => 0);
   New_Line;
end Sum_Series;

[edit] Aime

real
Invsqr(real n)
{
    return 1 / (n * n);
}
 
integer
main(void)
{
    integer i;
    real sum;
 
    sum = 0;
 
    i = 1;
    while (i < 1000) {
        sum += Invsqr(i);
        i += 1;
    }
 
    o_real(14, sum);
    o_byte('\n');
 
    return 0;
}

[edit] ALGOL 68

Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny
Works with: ELLA ALGOL 68 version Any (with appropriate job cards) - tested with release 1.8-8d
MODE RANGE = STRUCT(INT lwb, upb);
 
PROC sum = (PROC (INT)LONG REAL f, RANGE range)LONG REAL:(
  LONG REAL sum := LENG 0.0;
  FOR i FROM lwb OF range TO upb OF range DO
     sum := sum + f(i)
  OD;
  sum
);
 
test:(
  RANGE range = (1,100);
  PROC f = (INT x)LONG REAL: LENG REAL(1) / LENG REAL(x)**2;
  print(("Sum of f(x) from", lwb OF range, " to ",upb OF range," is ", SHORTEN sum(f,range),".", new line))
)

Output: 

Sum of f(x) from         +1 to        +100 is +1.63498390018489e  +0.

[edit] APL

      +/÷2*⍨⍳1000
1.64393

[edit] AutoHotkey

AutoHotkey allows the precision of floating point numbers generated by math operations to be adjusted via the SetFormat command. The default is 6 decimal places. 

SetFormat, FloatFast, 0.15
While A_Index <= 1000
 sum += 1/A_Index**2
MsgBox,% sum  ;1.643934566681554

[edit] AWK

$ awk 'BEGIN{for(i=1;i<=1000;i++)s+=1/(i*i);print s}'
1.64393

[edit] BASIC

Works with: QuickBasic version 4.5
FUNCTION s(x%)
   s = 1 / x ^ 2
END FUNCTION
 
FUNCTION sum(low%, high%)
   ret = 0
   FOR i = low TO high
      ret = ret + s(i)
   NEXT i
   sum = ret
END FUNCTION
PRINT sum(1, 1000)


[edit] BBC BASIC

      FOR i% = 1 TO 1000
        sum += 1/i%^2
      NEXT
      PRINT sum

[edit] Bracmat

( 0:?i
& 0:?S
& whl'(1+!i:~>1000:?i&!i^-2+!S:?S)
& out$!S
& out$(flt$(!S,10))
);

Output: 

8354593848314...../5082072010432.....  (1732 digits and a slash)
1,6439345667*10E0

[edit] Brat

p 1.to(1000).reduce 0 { sum, x | sum + 1.0 / x ^ 2 }  #Prints 1.6439345666816

[edit] C

#include <stdio.h>
 
double Invsqr(double n)
{
	return 1 / (n*n);
}
 
int main (int argc, char *argv[])
{
	int i, start = 1, end = 1000;
	double sum = 0.0;
 
	for( i = start; i <= end; i++)
		sum += Invsqr((double)i);           
 
	printf("%16.14f\n", sum);
 
	return 0;
}

[edit] C++

#include <iostream>
 
double f(double x);
 
int main()
{
    unsigned int start = 1;
    unsigned int end = 1000;
    double sum = 0;
 
    for( unsigned int x = start; x <= end; ++x )
    {
        sum += f(x);
    }
    std::cout << "Sum of f(x) from " << start << " to " << end << " is " << sum << std::endl;
    return 0;
}
 
 
double f(double x)
{
    return ( 1.0 / ( x * x ) );
}

[edit] C#

class Program
{
    static void Main(string[] args)
    {
        // Create and fill a list of number 1 to 1000
 
        List<double> myList = new List<double>();
        for (double i = 1; i < 1001; i++)
        {
            myList.Add(i);
        }
        // Calculate the sum of 1/x^2
 
        var sum = myList.Sum(x => 1/(x*x));
 
        Console.WriteLine(sum);
        Console.ReadLine();
    }
}

An alternative approach using Enumerable.Range() to generate the numbers. 

class Program
{
    static void Main(string[] args)
    {
        double sum = Enumerable.Range(1, 1000).Sum(x => 1.0 / (x * x));
 
        Console.WriteLine(sum);
        Console.ReadLine();
    }
}

[edit] CLIPS

(deffunction S (?x) (/ 1 (* ?x ?x)))
(deffunction partial-sum-S
  (?start ?stop)
  (bind ?sum 0)
  (loop-for-count (?i ?start ?stop) do
    (bind ?sum (+ ?sum (S ?i)))
  )
  (return ?sum)
)

Usage: 

CLIPS> (partial-sum-S 1 1000)
1.64393456668156

[edit] Clojure

(reduce + (map #(/ 1.0 % %) (range 1 1001)))

[edit] CoffeeScript

 
console.log [1..1000].reduce((acc, x) -> acc + (1.0 / (x*x)))
 

[edit] Common Lisp

(loop for x from 1 to 1000 summing (expt x -2))

[edit] D

[edit] More Procedural Style

import std.stdio, std.traits;
 
ReturnType!TF series(TF)(TF func, int end, int start=1) pure nothrow {
    typeof(return) sum = 0;
    foreach (i; start .. end + 1)
        sum += func(i);
    return sum;
}
 
void main() {
    writeln("Sum: ", series((int n) => 1.0L / (n ^^ 2), 1_000));
}
Output:
Sum: 1.64393

[edit] More functional Style

Same output. 

import std.stdio, std.algorithm, std.range;
 
auto series(alias F)(in int end, in int start=1) /*pure*/ {
    return iota(start, end + 1).map!F().reduce!q{a + b}();
}
 
void main() {
    writeln("Sum: ", series!q{1.0L / (a ^^ 2)}(1_000));
}

[edit] Delphi

 
 
unit Form_SumOfASeries_Unit;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;
 
type
  TFormSumOfASeries = class(TForm)
    M_Log: TMemo;
    B_Calc: TButton;
    procedure B_CalcClick(Sender: TObject);
  private
    { Private-Deklarationen }
  public
    { Public-Deklarationen }
  end;
 
var
  FormSumOfASeries: TFormSumOfASeries;
 
implementation
 
{$R *.dfm}
 
function Sum_Of_A_Series(_from,_to:int64):extended;
begin
  result:=0;
  while _from<=_to do
  begin
    result:=result+1.0/(_from*_from);
    inc(_from);
  end;
end;
 
procedure TFormSumOfASeries.B_CalcClick(Sender: TObject);
begin
  try
    M_Log.Lines.Add(FloatToStr(Sum_Of_A_Series(1,trunc(sqrt($7FFFFFFFFFFFFFFF-1)))));
  except
    M_Log.Lines.Add('Error');
  end;
end;
 
end.
 
 

Pi^2/6 = 1.644934066_848226436472415166646

result = 1,644934066_51266

[edit] DWScript

 
var s : Float;
for var i := 1 to 1000 do
   s += 1 / Sqr(i);
 
PrintLn(s);
 

[edit] E

pragma.enable("accumulator")
accum 0 for x in 1..1000 { _ + 1 / x ** 2 }

[edit] Eiffel

 
note
	description: "Compute the n-th term of a series"
 
class
	SUM_OF_SERIES_EXAMPLE
 
inherit
	MATH_CONST
 
create
	make
 
feature -- Initialization
 
	make
		local
			approximated, known: REAL_64
		do
			known := Pi^2 / 6
 
			approximated := sum_until (agent g, 1001)
			print ("%Nzeta function exact value: %N")
			print (known)
			print ("%Nzeta function approximated value: %N")
			print (approximated)
		end
 
feature -- Access
 
	g (k: INTEGER): REAL_64
			-- 'k'-th term of the serie
		require
			k_positive: k > 0
		do
			Result := 1 / (k * k)
		end
 
	sum_until (s: FUNCTION [ANY, TUPLE [INTEGER], REAL_64]; n: INTEGER): REAL_64
			-- sum of the 'n' first terms of 's'
		require
			n_positive: n > 0
			one_parameter: s.open_count = 1
		do
			Result := 0
			across 1 |..| n as it loop
				Result := Result + s.item ([it.item])
			end
		end
 
end
 
 

[edit] Erlang

lists:sum([1/math:pow(X,2) || X <- lists:seq(1,1000)]).

[edit] Euphoria

Works with: Euphoria version 4.0.0

This is based on the BASIC example. 

 
function s( atom x )
	return 1 / power( x, 2 )
end function 
 
function sum( atom low, atom high )
	atom ret = 0.0
	for i = low to high do
		ret = ret + s( i )
	end for
	return ret
end function
 
printf( 1, "%.15f\n", sum( 1, 1000 ) )

[edit] Factor

1000 [1,b] [ >float sq recip ] map-sum

[edit] Fantom

Within 'fansh': 

 
fansh> (1..1000).toList.reduce(0.0f) |Obj a, Int v -> Obj| { (Float)a + (1.0f/(v*v)) }
1.6439345666815615
 

[edit] Forth

: sum ( fn start count -- fsum )
  0e
  bounds do
    i s>d d>f dup execute f+
  loop drop ;
 
:noname ( x -- 1/x^2 ) fdup f* 1/f ;   ( xt )
1 1000 sum f.       \ 1.64393456668156
pi pi f* 6e f/ f.   \ 1.64493406684823

[edit] Fortran

In ISO Fortran 90 and later, use SUM intrinsic: 

real, dimension(1000) :: a = (/ (1.0/(i*i), i=1, 1000) /)
real :: result
 
result = sum(a);

[edit] F#

The following function will do the task specified. 

let rec f (x : float) = 
    match x with
        | 0. -> x
        | x -> (1. / (x * x)) + f (x - 1.)

In the interactive F# console, using the above gives: 

> f 1000. ;;
val it : float = 1.643934567

However this recursive function will run out of stack space eventually (try 100000). A tail-recursive implementation will not consume stack space and can therefore handle much larger ranges. Here is such a version: 

#light
let sum_series (max : float) =
    let rec f (a:float, x : float) = 
        match x with
            | 0. -> a
            | x -> f ((1. / (x * x) + a), x - 1.)
    f (0., max)
 
[<EntryPoint>]
let main args =
    let (b, max) = System.Double.TryParse(args.[0])
    printfn "%A" (sum_series max)
    0

This block can be compiled using fsc --target exe filename.fs or used interactively without the main function.

[edit] GAP

# We will compute the sum exactly
 
# Computing an approximation of a rationnal (giving a string)
# Value is truncated toward zero
Approx := function(x, d)
	local neg, a, b, n, m, s;
	if x < 0 then
		x := -x;
		neg := true;
	else
		neg := false;
	fi;
	a := NumeratorRat(x);
	b := DenominatorRat(x);
	n := QuoInt(a, b);
	a := RemInt(a, b);
	m := 10^d;
	s := "";
	if neg then
		Append(s, "-");
	fi;
	Append(s, String(n));
	n := Size(s) + 1;
	Append(s, String(m + QuoInt(a*m, b)));
	s[n] := '.';
	return s;
end;
 
a := Sum([1 .. 1000], n -> 1/n^2);;
Approx(a, 10);
"1.6439345666"
# and pi^2/6 is 1.6449340668, truncated to ten digits

[edit] Go

package main
 
import ("fmt"; "math")
 
func main() {
    fmt.Println("known:   ", math.Pi*math.Pi/6)
    sum := 0.
    for i := 1e3; i > 0; i-- {
        sum += 1 / (i * i)
    }
    fmt.Println("computed:", sum)
}

Output: 

known:    1.6449340668482264
computed: 1.6439345666815597

[edit] Groovy

Start with smallest terms first to minimize rounding error: 

println ((1000..1).collect { x -> 1/(x*x) }.sum())

Output: 

1.6439345654

[edit] Haskell

With a list comprehension: 

sum [1 / x ^ 2 | x <- [1..1000]]

With higher-order functions: 

sum $ map (\x -> 1 / x ^ 2) [1..1000]

In point-free style: 

(sum . map (1/) . map (^2)) [1..1000]

[edit] HicEst

REAL :: a(1000)
        a = 1 / $^2
        WRITE(ClipBoard, Format='F17.15') SUM(a) 
1.643934566681561

[edit] Icon and Unicon

procedure main()
   local i, sum
   sum := 0 & i := 0
   every sum +:= 1.0/((| i +:= 1 ) ^ 2) \1000
   write(sum)
end

or 

procedure main()
    every (sum := 0) +:= 1.0/((1 to 1000)^2)
    write(sum)
end

Note: The terse version requires some explanation. Icon expressions all return values or references if they succeed. As a result, it is possible to have expressions like these below: 

 
   x := y := 0   # := is right associative so, y is assigned 0, then x
   1 < x < 99    # comparison operators are left associative so, 1 < x returns x (if it is greater than 1), then x < 99 returns 99 if the comparison succeeds
   (sum := 0)    # returns a reference to sum which can in turn be used with augmented assignment +:=
 

[edit] IDL

print,total( 1/(1+findgen(1000))^2)

[edit] J

   NB. sum of reciprocals of squares of first thousand positive integers
   +/ % *: >: i. 1000
1.64393
 
   (*:o.1)%6       NB. pi squared over six, for comparison
1.64493
 
   1r6p2           NB.  As a constant (J has a rich constant notation)
1.64493

[edit] Java

public class Sum{
    public static double f(double x){
       return 1/(x*x);
    }
 
    public static void main(String[] args){
       double start = 1;
       double end = 1000;
       double sum = 0;
 
       for(double x = start;x <= end;x++) sum += f(x);
 
       System.out.println("Sum of f(x) from " + start + " to " + end +" is " + sum);
    }
}

[edit] JavaScript

function sum(a,b,fn) {
   var s = 0;
   for ( ; a <= b; a++) s += fn(a);
   return s;
}
 
 sum(1,1000, function(x) { return 1/(x*x) } )  // 1.64393456668156

[edit] Julia

julia> f(x) = sum(1/[1:x].^2)
# method added to generic function f
 
julia> f(1000)
1.6439345666815615

[edit] K

  ssr: +/1%_sqr
  ssr 1+!1000
1.643935

[edit] Lang5

1000 iota 1 + 1 swap / 2 ** '+ reduce .

[edit] Liberty BASIC

 
for i =1 to 1000
  sum =sum +1 /( i^2)
next i
 
print sum
 
end
 

[edit] 

to series :fn :a :b
  localmake "sigma 0
  for [i :a :b] [make "sigma :sigma + invoke :fn :i]
  output :sigma
end
to zeta.2 :x
  output 1 / (:x * :x)
end
print series "zeta.2 1 1000
make "pi (radarctan 0 1) * 2
print :pi * :pi / 6

[edit] Lua

 
sum = 0
for i = 1, 1000 do sum = sum + 1/i^2 end
print(sum)
 

[edit] Lucid

series = ssum asa  n >= 1000
   where
         num = 1 fby num + 1;
         ssum = ssum + 1/(num * num)
   end;

[edit] Mathematica

This is the straightforward solution of the task: 

Sum[1/x^2, {x, 1, 1000}]

However this returns a quotient of two huge integers (namely the exact sum); to get a floating point approximation, use N: 

N[Sum[1/x^2, {x, 1, 1000}]]

or better: 

NSum[1/x^2, {x, 1, 1000}]

Which gives a higher or equal accuracy/precision. Alternatively, get Mathematica to do the whole calculation in floating point by using a floating point value in the formula: 

Sum[1./x^2, {x, 1, 1000}]

Other ways include (exact, approximate,exact,approximate): 

Total[Table[1/x^2, {x, 1, 1000}]]   
Total[Table[1./x^2, {x, 1, 1000}]]
Plus@@Table[1/x^2, {x, 1, 1000}]
Plus@@Table[1./x^2, {x, 1, 1000}]

[edit] MATLAB

   sum([1:1000].^(-2))

[edit] Maxima

(%i45) sum(1/x^2, x, 1, 1000);
       835459384831496894781878542648[806 digits]396236858699094240207812766449
(%o45) ------------------------------------------------------------------------
       508207201043258126178352922730[806 digits]886537101453118476390400000000
 
(%i46) sum(1/x^2, x, 1, 1000),numer;
(%o46) 1.643934566681561

[edit] MAXScript

total = 0
for i in 1 to 1000 do
(
    total += 1.0 / pow i 2
)
print total

[edit] МК-61/52

0	П0	П1	ИП1	1	+	П1	x^2	1/x	ИП0
+	П0	ИП1	1	0	0	0	-	x>=0	03
ИП0	С/П

[edit] MMIX

x	IS	$1	% flt calculations
y	IS	$2	%   id
z	IS	$3	% z = sum series
t	IS	$4	% temp var
 
	LOC	Data_Segment
	GREG	@
BUF	OCTA	0,0,0		% print buffer
 
	LOC	#1000
	GREG	@

// print floating point number in scientific format: 0.xxx...ey.. 
// most of this routine is adopted from:
// http://www.pspu.ru/personal/eremin/emmi/rom_subs/printreal.html
// float number in z
	GREG	@
NaN	BYTE	"NaN..",0
NewLn	BYTE	#a,0
1H	LDA	x,NaN
	TRAP	0,Fputs,StdOut
	GO	$127,$127,0
 
prtFlt	FUN	x,z,z		% test if z == NaN
	BNZ	x,1B
	CMP	$73,z,0		% if necessary remember it is neg
	BNN	$73,4F
Sign	BYTE	'-'
	LDA	$255,Sign
	TRAP	0,Fputs,StdOut
	ANDNH	z,#8000		% make number pos
// normalizing float number
4H	SETH	$74,#4024	% initialize mulfactor = 10.0
	SETH	$73,#0023
	INCMH	$73,#86f2
	INCML	$73,#6fc1	%
	FLOT	$73,$73		% $73 = float 10^16
	SET	$75,16		% set # decimals to 16
8H	FCMP	$72,z,$73	% while z >= 10^16 do
	BN	$72,9F		% 
	FDIV	z,z,$74		%  z = z / 10.0 
	ADD	$75,$75,1	%  incr exponent
	JMP	8B		% wend
9H	FDIV	$73,$73,$74	% 10^16 / 10.0
5H	FCMP	$72,z,$73	% while z < 10^15 do 
	BNN	$72,6F
	FMUL	z,z,$74		%  z = z * 10.0
	SUB	$75,$75,1	%  exp = exp - 1
	JMP	5B
NulPnt	BYTE	'0','.',#00
6H	LDA	$255,NulPnt	% print '0.' to StdOut
	TRAP	0,Fputs,StdOut
	FIX	z,0,z		% convert float z to integer 
// print mantissa
0H	GREG	#3030303030303030
	STO	0B,BUF
	STO	0B,BUF+8	% store print mask in buffer
	LDA	$255,BUF+16	% points after LSD
				% repeat
2H	SUB	$255,$255,1	%   move pointer down
	DIV	z,z,10		%   (q,r) = divmod z 10 
	GET	t,rR		%   get remainder
	INCL	t,'0'		%   convert to ascii digit
	STBU	t,$255,0	%   store digit in buffer
	BNZ	z,2B		% until q == 0
	TRAP	0,Fputs,StdOut	% print mantissa
Exp	BYTE	'e',#00
	LDA	$255,Exp	% print 'exponent' indicator
	TRAP	0,Fputs,StdOut
// print exponent
0H	GREG	#3030300000000000
	STO	0B,BUF
	LDA	$255,BUF+2	% store print mask in buffer
	CMP	$73,$75,0	% if exp neg then place - in buffer
	BNN	$73,2F
ExpSign	BYTE	'-'	
	LDA	$255,ExpSign
	TRAP	0,Fputs,StdOut
	NEG	$75,$75		% make exp positive
2H	LDA	$255,BUF+3	% points after LSD
				% repeat
3H	SUB	$255,$255,1	%   move pointer down
	DIV	$75,$75,10	%   (q,r) = divmod exp 10
	GET	t,rR
	INCL	t,'0'
	STBU	t,$255,0	%   store exp. digit in buffer
	BNZ	$75,3B		% until q == 0
	TRAP	0,Fputs,StdOut	% print exponent
	LDA	$255,NewLn
	TRAP	0,Fputs,StdOut	% do a NL
	GO	$127,$127,0	% return
 
i  IS $5 ;iu IS $6
Main	SET	iu,1000
	SETH	y,#3ff0     y = 1.0
	SETH	z,#0000     z = 0.0
	SET	i,1          for (i=1;i<=1000; i++ ) {
1H	FLOT	x,i           x = int i
	FMUL	x,x,x         x = x^2
	FDIV	x,y,x         x = 1 / x
	FADD	z,z,x         s = s + x
	ADD	i,i,1
	CMP	t,i,iu
	PBNP	t,1B         } z = sum
	GO	$127,prtFlt  print sum --> StdOut
	TRAP	0,Halt,0

Output: 

~/MIX/MMIX/Rosetta> mmix sumseries
0.1643934566681562e1

[edit] Modula-3

Modula-3 uses D0 after a floating point number as a literal for LONGREAL. 

MODULE Sum EXPORTS Main;
 
IMPORT IO, Fmt, Math;
 
VAR sum: LONGREAL := 0.0D0;
 
PROCEDURE F(x: LONGREAL): LONGREAL =
  BEGIN
    RETURN 1.0D0 / Math.pow(x, 2.0D0);
  END F;
 
BEGIN
  FOR i := 1 TO 1000 DO
    sum := sum + F(FLOAT(i, LONGREAL));
  END;
  IO.Put("Sum of F(x) from 1 to 1000 is ");
  IO.Put(Fmt.LongReal(sum));
  IO.Put("\n");
END Sum.

Output: 

Sum of F(x) from 1 to 1000 is 1.6439345666815612

[edit] MUMPS

 
SOAS(N)
 NEW SUM,I SET SUM=0
 FOR I=1:1:N DO
 .SET SUM=SUM+(1/((I*I)))
 QUIT SUM
 

This is an extrinsic function so the usage is: 

USER>SET X=$$SOAS^ROSETTA(1000) WRITE X
1.643934566681559806

[edit] Nial

|sum (1 / power (count 1000) 2)
=1.64393

[edit] Objeck

 
bundle Default {
  class SumSeries {
    function : Main(args : String[]) ~ Nil {
      DoSumSeries();
    }
 
    function : native : DoSumSeries() ~ Nil {
      start := 1;
      end := 1000;
 
      sum := 0.0;
 
      for(x : Float := start; x <= end; x += 1;) {
        sum += f(x);
      };
 
      IO.Console->GetInstance()->Print("Sum of f(x) from ")->Print(start)->Print(" to ")->Print(end)->Print(" is ")->PrintLine(sum);
    }
 
    function : native : f(x : Float) ~ Float {
      return 1.0 / (x * x);
    }
  }
}
 

[edit] OCaml

let sum a b fn =
  let result = ref 0. in
  for i = a to b do
    result := !result +. fn i
  done;
  !result
# sum 1 1000 (fun x -> 1. /. (float x ** 2.))
- : float = 1.64393456668156124

or in a functional programming style: 

let sum a b fn =
  let rec aux i r =
    if i > b then r
    else aux (succ i) (r +. fn i)
  in
  aux a 0.
;;

Simple recursive solution: 

let rec sum n = if n < 1 then 0.0 else sum (n-1) +. 1.0 /. float (n*n)
in sum 1000

[edit] Octave

Given a vector, the sum of all its elements is simply sum(vector); a range can be generated through the range notation: sum(1:1000) computes the sum of all numbers from 1 to 1000. To compute the requested series, we can simply write: 

sum(1 ./ [1:1000] .^ 2)


[edit] OpenEdge/Progress

Conventionally like elsewhere: 

DEF VAR dcResult AS DECIMAL NO-UNDO.
DEF VAR n AS INT NO-UNDO.
 
DO n = 1 TO 1000 :
  dcResult = dcResult + 1 / (n * n)  .
END.
 
DISPLAY dcResult .

or like this: 

DEF VAR n AS INT NO-UNDO.
 
REPEAT n = 1 TO 1000 :
  ACCUMULATE 1 / (n * n) (total).
END.
 
DISPLAY ( ACCUM total 1 / (n * n) )  .

[edit] Oz

With higher-order functions: 

declare
  fun {SumSeries S N}
     {FoldL {Map {List.number 1 N 1} S}
      Number.'+' 0.}
  end
 
  fun {S X}
     1. / {Int.toFloat X*X}
  end
in
  {Show {SumSeries S 1000}}

Iterative: 

  fun {SumSeries S N}
     R = {NewCell 0.}
  in
     for I in 1..N do
        R := @R + {S I}
     end
     @R
  end

[edit] PARI/GP

Exact rational solution: 

sum(n=1,1000,1/n^2)

Real number solution (accurate to 3\cdot10^{-36} at standard precision): 

sum(n=1,1000,1./n^2)

Approximate solution (accurate to 9\cdot10^{-11} at standard precision): 

zeta(2)-intnum(x=1000.5,[1],1/x^2)

or 

zeta(2)-1/1000.5

[edit] Pascal

Program SumSeries;
 
var
  S: double;
  i: integer;
 
function f(number: integer): double;
begin
  f := 1/(number*number);
end;
 
begin
  S := 0;
  for i := 1 to 1000 do
    S := S + f(i);
  writeln('The sum of 1/x^2 from 1 to 1000 is: ', S:10:8);
  writeln('Whereas pi^2/6 is:                  ', pi*pi/6:10:8);
end.

Output: 

The sum of a series from 1 to 1000 is: 1.64393457
Whereas pi^2/6 is:                     1.64493407

[edit] Perl

my $sum = 0;
$sum += 1 / $_ ** 2 foreach 1..1000;
print "$sum\n";

or 

use List::Util qw(reduce);
$sum = reduce { $a + 1 / $b ** 2 } 0, 1..1000;
print "$sum\n";

An other way of doing it is to define the series as a closure: 

my $S = do { my ($sum, $k); sub { $sum += 1/++$k**2 } };
my @S = map &$S, 1 .. 1000;
print $S[-1];

[edit] Perl 6

(Some of these work with rakudo, and others with niecza. Eventually they'll all work everywhere...)

In general, the $nth partial sum of a series whose terms are given by a unary function &f is 

[+] map &f, 1 .. $n

So what's needed in this case is 

say [+] map { 1 / $^n**2 }, 1 .. 1000;

Or, using the "hyper" metaoperator to vectorize, we can use a more "point free" style while keeping traditional precedence: 

say [+] 1 «/« (1..1000) »**» 2;

Or we can use the X "cross" metaoperator, which is convenient even if one side or the other is a scalar. In this case, we demonstrate a scalar on either side: 

say [+] 1 X/ (1..1000 X** 2);

Note that cross ops are parsed as list infix precedence rather than using the precedence of the base op as hypers do. Hence the difference in parenthesization.

In a lazy language like Perl 6, it's generally considered a stronger abstraction to write the correct infinite sequence, and then take the part of it you're interested in. Here we define an infinite sequence of partial sums (by adding a backslash into the reduction to make it look "triangular"), then take the 1000th term of that: 

constant @x = [\+] 0, { 1 / ++(state $n) ** 2 } ... *;
say @x[1000];  # prints 1.64393456668156

Note that infinite constant sequences can be lazily generated in Perl 6, or this wouldn't work so well...

A cleaner style is to combine these approaches with a more FP look: 

constant ζish = [\+] map -> \𝑖 { 1 / 𝑖**2 }, 1..*;
say ζish[1000];

Perhaps the cleanest way is to just define the zeta function and evaluate it for s=2: 

sub ζ($s) { [\+] map -> \n { 1 / n**$s }, 1..* }
say ζ(2)[1000];

Notice how the thus-defined zeta function returns a lazy list of approximated values, which is arguably the closest we can get from the mathematical definition.

Finally, if list comprehensions are your hammer, you can nail it this way: 

say [+] (1 / $_**2 for 1..1000);

That's fine for a single result, but if you're going to be evaluating the sequence multiple times, you don't want to be recalculating the sum each time, so it's more efficient to define the sequence as a constant to let the run-time automatically cache those values already calculated.

[edit] PicoLisp

(scl 9)  # Calculate with 9 digits precision
 
(let S 0
   (for I 1000
      (inc 'S (*/ 1.0 (* I I))) )
   (prinl (round S 6)) )  # Round result to 6 digits

Output: 

1.643935

[edit] Pike

array(int) x = enumerate(1000,1,1);
`+(@(1.0/pow(x[*],2)[*]));
Result: 1.64393

[edit] PHP

<?php
/**
 * @author Elad Yosifon <[email protected]>
 */
 
/**
 * @param int $n
 * @param int $k
 * @return float|int
 */
function sum_of_a_series($n,$k)
{
	$sum_of_a_series = 0;
	for($i=$k;$i<=$n;$i++)
	{
		$sum_of_a_series += (1/($i*$i));
	}
	return $sum_of_a_series;
}
 
echo sum_of_a_series(1000,1);
 
Output:
1.6439345666816

[edit] PL/I

 
/* sum the first 1000 terms of the series 1/n**2. */
s = 0;
 
do i = 1000 to 1 by -1;
   s = s + 1/float(i**2);
end;
 
put skip list (s);
 

Output: 

1.64393456668155980E+0000

[edit] Pop11

lvars s = 0, j;
for j from 1 to 1000 do
    s + 1.0/(j*j) -> s;
endfor;
 
s =>

[edit] PostScript

 
/aproxriemann{
/x exch def
/i 1 def
/sum 0 def
x{
/sum sum i -2 exp add def
/i i 1 add def
}repeat
sum ==
}def
 
1000 aproxriemann
 

Output: 

 
1.64393485
 
Library: initlib
 
% using map
[1 1000] 1 range {dup * 1 exch div} map 0 {+} fold
 
% just using fold
[1 1000] 1 range 0 {dup * 1 exch div +}  fold
 

[edit] PowerShell

$x = 1..1000 `
       | ForEach-Object { 1 / ($_ * $_) } `
       | Measure-Object -Sum
Write-Host Sum = $x.Sum

[edit] Prolog

Works with SWI-Prolog. 

sum(S) :-
        findall(L, (between(1,1000,N),L is 1/N^2), Ls),
        sumlist(Ls, S).
 

Ouptput : 

?- sum(S).
S = 1.643934566681562.

[edit] PureBasic

Define i, sum.d
 
For i=1 To 1000
  sum+1.0/(i*i)
Next i
 
Debug sum

Answer = 1.6439345666815615

[edit] Python

print ( sum(1.0 / (x * x) for x in range(1, 1001)) )

[edit] R

print( sum( 1/seq(1000)^2 ) )

[edit] Racket

A solution using Typed Racket: 

 
#lang typed/racket
 
(: S : Natural -> Real)
(define (S n)
  (for/sum: : Real ([k : Natural (in-range 1 (+ n 1))])
    (/ 1.0 (* k k))))
 

[edit] Raven

0 1 1000 1 range each 1.0 swap dup * / +
"%g\n" print
Output:
1.64393

Raven uses a 32 bit float, so precision limits the accuracy of the result for large iterations.

[edit] REXX

[edit] sums specific terms

/*REXX program sums the first   N   terms of    1/(i**2),    i=1 ──►  N.*/
parse arg N D .                        /*maybe get num of terms, digits.*/
if N=='' | N==','  then N=1000         /*Not specified?  Use the default*/
if D==''           then D=60           /* "      "        "   "     "   */
numeric digits D                       /*use D digits: 9 is default for */
                                       /*REXX, 60 is this pgm's default.*/
w = length(N)                          /*use max  width for nice output.*/
sum = 0                                /*initialize the sum to zero.    */
          do j=1  for N                /*compute for   N   terms.       */
          sum = sum   +   1 / j**2     /*add another term to the sum.   */
          end   /*j*/
 
say  'The sum of'   right(N,w)   "terms is:"   sum
                                       /*stick a fork in it, we're done.*/

output when using the default input 

The sum of 1000 terms is: 1.64393456668155980313905802382221558965210344649368531671713

[edit] sums with running total

/*REXX program sums the first   N   terms of    1/(i**2),    i=1 ──►  N.*/
parse arg N D .                        /*maybe get num of terms, digits.*/
if N=='' | N==','  then N=1000         /*Not specified?  Use the default*/
if D==''           then D=60           /* "      "        "   "     "   */
numeric digits D                       /*use D digits: 9 is default for */
                                       /*REXX, 60 is this pgm's default.*/
w = length(N)                          /*use max  width for nice output.*/
sum=0                                  /*initialize the sum to zero.    */
      do j=1  for N                    /*compute for   N   terms.       */
      sum = sum   +   1 / j**2         /*add another term to the sum.   */
      if   left(j,1)\==1  then iterate /*does J start with a one ?      */
      if  right(j,1)\==0  then iterate /*  "  "  end   "   " zero ?     */
      if substr(j,2)\= 0  then iterate /*  "  "   "    "   all zeroes ? */
      say  'The sum of'   right(j,w)   "terms is:"   sum   /*display it.*/
      end   /*j*/
                                       /*stick a fork in it, we're done.*/

output when using the input of: 1000000000 

The sum of         10 terms is: 1.54976773116654069035021415973796926177878558830939783320736
The sum of        100 terms is: 1.63498390018489286507716949818032376668332170003126381385307
The sum of       1000 terms is: 1.64393456668155980313905802382221558965210344649368531671713
The sum of      10000 terms is: 1.64483407184805976980608183331031090353799751949684175308996
The sum of     100000 terms is: 1.64492406689822626980574850331269185564752132981156034248806
The sum of    1000000 terms is: 1.64493306684872643630574849997939185588561654406394129491321
The sum of   10000000 terms is: 1.64493396684823143647224849997935852288561656787346272343397
The sum of  100000000 terms is: 1.64493405684822648647241499997935852255228656787346510441026
The sum of 1000000000 terms is: 1.64493406584822643697241516647935852255228323457346510444171

output from a calculator {π²/6, using 60 digits) showing the correct number of digits [superscript was added after-the-fact]: 

1.64493406684822643647241516664602518921894990120679843773556

[edit] sums with running significance

This is a technique to show a running significance (based on the previous calculation).

If the old REXX variable would be set to 1.64 (instead of 1), the first noise digits could be
bypassed to make the display cleaner. 

/*REXX program sums the first   N   terms of    1/(i**2),    i=1 ──►  N.*/
parse arg N D .                        /*optional num of terms,  digits.*/
if N=='' | N==','  then N=1000         /*Not specified?  Use the default*/
if D==''           then D=60           /* "      "        "   "     "   */
@sig = 'The significant sum of'        /*literal used in  SAY statement.*/
numeric digits D                       /*use D digits: 9 is default for */
                                       /*REXX, 60 is this pgm's default.*/
w = length(N)                          /*use max  width for nice output.*/
sum = 0                                /*initialize the  SUM  to zero.  */
old = 1                                /*the SUM to compared to the NEW.*/
p = 0                                  /*significant precision so far.  */
      do j=1  for N                    /*compute for   n   terms.       */
      sum = sum   +   1 / j**2         /*add another term to the sum.   */
      c = compare(sum,old)             /*see how we're doing with prec. */
      if c>p  then do                  /*Got another significant digit? */
                   say  @sig   right(j,w)    "terms is:"    left(sum,c)
                   p = c               /*the new significant precision. */
                   end
      old = sum                        /*use "old" sum for next compare.*/
      end   /*j*/
say
say  'The sum of'   right(N,w)   "terms is:"   /*display sum's preamble.*/
say  sum                               /*display the sum on its own line*/
                                       /*stick a fork in it, we're done.*/

output when using the input of: 35000000 100 

The significant sum of        2 terms is: 1.
The significant sum of        3 terms is: 1.3
The significant sum of        5 terms is: 1.46
The significant sum of       14 terms is: 1.575
The significant sum of       34 terms is: 1.6159
The significant sum of      110 terms is: 1.63588
The significant sum of      328 terms is: 1.641889
The significant sum of     1024 terms is: 1.6439579
The significant sum of     3207 terms is: 1.64462229
The significant sum of    10043 terms is: 1.644834499
The significant sum of    31782 terms is: 1.6449026029
The significant sum of   100314 terms is: 1.64492409819
The significant sum of   316728 terms is: 1.644930909569
The significant sum of  1000853 terms is: 1.6449330677009
The significant sum of  3163463 terms is: 1.64493375073899
The significant sum of 10001199 terms is: 1.644933966860219
The significant sum of 31627592 terms is: 1.6449340352302649

The sum of 35000000 terms is:
1.644934038276798273207105156927852205740478629316117966926591883437164764834567731984252290795163298

One can see a pattern in the number of significant digits computed based on the number of terms used.
(See a discussion in the talk section.)

[edit] RLaB

 
>> sum( (1 ./ [1:1000]) .^ 2 ) - const.pi^2/6
-0.000999500167
 

[edit] Ruby

puts (1..1000).inject(0) {|sum, x| sum + 1.0 / x ** 2}

[edit] Run BASIC

 
for i =1 to 1000
  sum = sum + 1 /( i^2)
next i
print sum

[edit] SAS

data _null_;
s=0;
do n=1 to 1000;
   s+1/n**2;        /* s+x is synonym of s=s+x */
end;
e=s-constant('pi')**2/6;
put s e;
run;

[edit] Scala

scala> 1 to 1000 map (x => 1.0 / (x * x)) sum
res30: Double = 1.6439345666815615

[edit] Scheme

(define (sum a b fn)
  (do ((i a (+ i 1))
       (result 0 (+ result (fn i))))
      ((> i b) result)))
 
(sum 1 1000 (lambda (x) (/ 1 (* x x)))) ; fraction
(exact->inexact (sum 1 1000 (lambda (x) (/ 1 (* x x))))) ; decimal

More idiomatic way (or so they say) by tail recursion: 

(define (invsq f to)
  (let loop ((f f) (s 0))
    (if (> f to)
      s
      (loop (+ 1 f) (+ s (/ 1 f f))))))
 
;; whether you get a rational or a float depends on implementation
(invsq 1 1000) ; 835459384831...766449/50820...90400000000
(exact->inexact (invsq 1 1000)) ; 1.64393456668156

[edit] Seed7

$ include "seed7_05.s7i";
  include "float.s7i";
 
const func float: invsqr (in float: n) is
  return 1.0 / n**2;
 
const proc: main is func
  local
    var integer: i is 0;
    var float: sum is 0.0;
  begin
    for i range 1 to 1000 do
      sum +:= invsqr(flt(i));
    end for;
    writeln(sum digits 6 lpad 8);
  end func;

[edit] Slate

Manually coerce it to a float, otherwise you will get an exact (and slow) answer: 

((1 to: 1000) reduce: [|:x :y | x + (y squared reciprocal as: Float)]).

[edit] Smalltalk

( (1 to: 1000) fold: [:sum :aNumber |
  sum + (aNumber squared reciprocal) ] ) asFloat displayNl.

[edit] SQL

CREATE TABLE SEQUENCE (
	n REAL
)
 
INSERT INTO SEQUENCE VALUES (0)
INSERT INTO SEQUENCE VALUES (1)
INSERT INTO SEQUENCE SELECT 2+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 4+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 8+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 16+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 32+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 64+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 128+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 256+n FROM SEQUENCE
INSERT INTO SEQUENCE SELECT 512+n FROM SEQUENCE
 
SELECT SUM(1/n) FROM SEQUENCE WHERE n>=1 AND n<=1000

This produces a result set with a value which is something like 7.48547092382796

[edit] Standard ML

 
(* 1.64393456668 *)
List.foldl op+ 0.0 (List.tabulate(1000, fn x => 1.0 / Math.pow(real(x + 1),2.0)))
 

[edit] Tcl

Works with: Tcl version 8.5
package require Tcl 8.5
 
proc partial_sum {func - start - stop} {
    for {set x $start; set sum 0} {$x <= $stop} {incr x} {
        set sum [expr {$sum + [apply $func $x]}]
    }
    return $sum
}
 
set S {x {expr {1.0 / $x**2}}}
 
partial_sum $S from 1 to 1000 ;# => 1.6439345666815615
Library: Tcllib (Package: struct::list)
package require Tcl 8.5
package require struct::list
 
proc sum_of {lambda nums} {
    struct::list fold [struct::list map $nums [list apply $lambda]] 0 ::tcl::mathop::+
}
 
sum_of $S [range 1 1001] ;# ==> 1.6439345666815615

The helper range procedure is: 

# a range command akin to Python's
proc range args {
    foreach {start stop step} [switch -exact -- [llength $args] {
        1 {concat 0 $args 1}
        2 {concat   $args 1}
        3 {concat   $args  }
        default {error {wrong # of args: should be "range ?start? stop ?step?"}}
    }] break
    if {$step == 0} {error "cannot create a range when step == 0"}
    set range [list]
    while {$step > 0 ? $start < $stop : $stop < $start} {
        lappend range $start
        incr start $step
    }
    return $range
}

[edit] TI-89 BASIC

∑(1/x^2,x,1,1000)

[edit] UnixPipes

term() {
   b=$1;res=$2
   echo "scale=5;1/($res*$res)+$b" | bc
}
 
sum() {
  (read B; res=$1;
  test -n "$B" && (term $B $res) || (term 0 $res))
}
 
fold() {
  func=$1
  (while read a ; do
      fold $func | $func $a
  done)
}
 
(echo 3; echo 1; echo 4) | fold sum

[edit] Unicon

See Icon.

[edit] Ursala

The expression plus:-0. represents a function returning the sum of any given list of floating point numbers, or zero if it's empty, using the built in reduction operator, :-, and the binary addition function, plus. The rest the expression constructs the series by inverting the square of each number in the list from 1 to 1000. 

#import flo
#import nat
 
#cast %e
 
total = plus:-0 div/*1. sqr* float*t iota 1001

output: 

1.643935e+00

[edit] Vala

 
public static void main(){
	int i, start = 1, end = 1000;
	double sum = 0.0;
 
	for(i = start; i<= end; i++)
		sum += (1 / (double)(i * i));
 
	stdout.printf("%s\n", sum.to_string());
}
 

Output: 

1.6439345666815615

[edit] XPL0

code CrLf=9;  code real RlOut=48;
int  X;  real S;
[S:= 0.0;
for X:= 1 to 1000 do S:= S + 1.0/float(X*X);
RlOut(0, S);  CrLf(0);
]

Output: 

    1.64393

[edit] Yorick

(1./indgen(1:1000)^2)(sum)
Retrieved from "http://rosettacode.org/mw/index.php?title=Sum_of_a_series&oldid=159737"
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox