Continued fraction: Difference between revisions

 

=={{header|11l}}==

V r = 0.0

L n > 0

print(calc(n -> I n > 0 {2} E 1, n -> 1))

print(calc(n -> I n > 0 {n} E 2, n -> I n > 1 {n - 1} E 1))

 

=={{header|Ada}}==

 

Generic function for estimating continued fractions:

type Scalar is digits <>;

 

with function A (N : in Natural) return Natural;

with function B (N : in Positive) return Natural;

 

function A (N : in Natural) return Scalar is (Scalar (Natural'(A (N))));

function B (N : in Positive) return Scalar is (Scalar (Natural'(B (N))));

end loop;

return A (0) + Fraction;

Test program using the function above to estimate the square root of 2, Napiers constant and pi:

 

with Continued_Fraction;

Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;

Put (Pi.Estimate (10000), Exp => 0); New_Line;

===Using only Ada 95 features===

This example is exactly the same as the preceding one, but implemented using only Ada 95 features.

type Scalar is digits <>;

 

with function A (N : in Natural) return Natural;

with function B (N : in Positive) return Natural;

 

function A (N : in Natural) return Scalar is

begin

end loop;

return A (0) + Fraction;

 

with Continued_Fraction_Ada95;

Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;

Put (Pi.Estimate (10000), Exp => 0); New_Line;

{{out}}

<pre> 1.41421356237310

=={{header|ALGOL 68}}==

{{works with|Algol 68 Genie|2.8}}

PROC cf = (INT steps, PROC (INT) INT a, PROC (INT) INT b) REAL:

BEGIN

print (("Napier: ", cf(precision, anap, bnap), newline));

print (("Pi: ", cf(precision, api, bpi)))

{{out}}

<pre>

Pi: +3.14159265358954e +0

</pre>

 

=={{header|ATS}}==

A fairly direct translation of the [[#C|C version]] without using advanced features of the type system:

"share/atspre_staload.hats"

//

println! ("napier = ", calc(napier, 100));

println! (" pi = ", calc( pi , 100));

 

=={{header|AutoHotkey}}==

{

return n?2.0:1.0

}

 

Output with Autohotkey v1 (currently 1.1.16.05):

e = 2.718282

Output with Autohotkey v2 (currently alpha 56):

e = 2.7182818284590455

Note the far superiour accuracy of v2.

 

=={{header|Axiom}}==

Axiom provides a ContinuedFraction domain:

get continuedFraction(1, repeating [1], repeating [2]) :: Float

get continuedFraction(2, cons(1,[i for i in 1..]), [i for i in 1..]) :: Float

Type: Float

 

 

(3) 3.1415926538 39792926

The value for <math>\pi</math> has an accuracy to only 9 decimal places after 1000 iterations, with an accuracy to 12 decimal places after 10000 iterations.

 

We could re-implement this, with the same output:

n=1 => initial

temp := 0

cf(1, repeating [1], repeating [2], 1000) :: Float

cf(2, cons(1,[i for i in 1..]), [i for i in 1..], 1000) :: Float

 

=={{header|BBC BASIC}}==

{{works with|BBC BASIC for Windows}}

@% = &1001010

expr$ += STR$(EVAL(a$)) + "+" + STR$(EVAL(b$)) + "/("

UNTIL LEN(expr$) > (65500 - N)

{{out}}

<pre>

=={{header|C}}==

{{works with|ANSI C}}

#include <stdio.h>

 

 

return 0;

{{out}}

<pre> 1.414213562

=={{header|C sharp|C#}}==

{{trans|Java}}

using System.Collections.Generic;

 

}

}

{{out}}

<pre>1.4142135623731

 

=={{header|C++}}==

#include <iostream>

#include <tuple>

<< calc(pi, 10000) << '\n'

<< std::setprecision(old_prec); // reset precision

{{out}}

<pre>

 

Functions don't take other functions as arguments, so I wrapped them in a dummy record each.

var r = 0.0;

 

writeln(calc(new Sqrt2(), n));

writeln(calc(new Napier(), n));

 

=={{header|Clojure}}==

(defn cfrac

[a b n]

(def e (cfrac #(if (zero? %) 2.0 %) #(if (= 1 %) 1.0 (double (dec %))) 100))

(def pi (cfrac #(if (zero? %) 3.0 6.0) #(let [x (- (* 2.0 %) 1.0)] (* x x)) 900000))

{{out}}

<pre>

{{works with|GnuCOBOL}}

 

program-id. show-continued-fractions.

 

goback.

end program pi-beta.

 

{{out}}

 

=={{header|CoffeeScript}}==

# a0 + b1 / (a1 + b2 / (a2 + b3 / ...

continuous_fraction = (f) ->

x = 2*n - 1

x * x

{{out}}

<pre>

=={{header|Common Lisp}}==

{{trans|C++}}

(let ((temp 0))

(loop for n1 from n downto 1

(* (1- (* 2 n))

(1- (* 2 n))))) 10000)

{{out}}

<pre>sqrt(2) = 1.4142135623730947d0

 

=={{header|D}}==

 

FP calc(FP, F)(in F fun, in int n) pure nothrow if (isCallable!F) {

calc!real(&fNapier, 200).print;

calc!real(&fPi, 200).print;

{{out}}

<pre>1.4142135623730950487

2.7182818284590452354

3.1415926228048469486</pre>

 

=={{header|Erlang}}==

-module(continued).

-compile([export_all]).

test_nappier (N) ->

continued_fraction(fun nappier_a/1,fun nappier_b/1,N).

 

{{out}}

29> continued:test_pi(1000).

3.141592653340542

31> continued:test_nappier(1000).

2.7182818284590455

 

=={{header|F_Sharp|F#}}==

===The Functions===

// I provide four functions:-

// cf2S general purpose continued fraction to sequence of float approximations

let cfSqRt n=(cf2S (fun()->n-1M) (let mutable n=false in fun()->if n then 2M else (n<-true; 1M)))

let π n=let mutable π=n in (fun ()->match π with h::t->π<-t; h |_->9999999999999999999999999999M)

===The Tasks===

cfSqRt 2M |> Seq.take 10 |> Seq.pairwise |> Seq.iter(fun(n,g)->printfn "%1.14f < √2 < %1.14f" (min n g) (max n g))

{{out}}

<pre>

1.41421355164605 < √2 < 1.41421362489487

</pre>

cfSqRt 0.25M |> Seq.take 30 |> Seq.iter (printfn "%1.14f")

{{out}}

<pre>

0.50000000000000

</pre>

let aπ()=let mutable n=0M in (fun ()->n<-n+1M;let b=n+n-1M in b*b)

let bπ()=let mutable n=true in (fun ()->match n with true->n<-false;3M |_->6M)

cf2S (aπ()) (bπ()) |> Seq.take 10 |> Seq.pairwise |> Seq.iter(fun(n,g)->printfn "%1.14f < π < %1.14f" (min n g) (max n g))

{{out}}

<pre>

3.14140671849650 < π < 3.14183961892940

</pre>

let pi = π [3M;7M;15M;1M;292M;1M;1M;1M;2M;1M;3M;1M;14M;2M;1M;1M;2M;2M;2M;2M]

cN2S pi |> Seq.take 10 |> Seq.pairwise |> Seq.iter(fun(n,g)->printfn "%1.14f < π < %1.14f" (min n g) (max n g))

{{out}}

<pre>

3.14159265358939 < π < 3.14159265359140

</pre>

let ae()=let mutable n=0.5M in (fun ()->match n with 0.5M->n<-0M; 1M |_->n<-n+1M; n)

let be()=let mutable n=0.5M in (fun ()->match n with 0.5M->n<-0M; 2M |_->n<-n+1M; n)

cf2S (ae()) (be()) |> Seq.take 10 |> Seq.pairwise |> Seq.iter(fun(n,g)->printfn "%1.14f < e < %1.14f" (min n g) (max n g))

{{out}}

<pre>

2.71828165766640 < e < 2.71828184277783

2.71828182735187 < e < 2.71828184277783

</pre>

 

=={{header|Factor}}==

''cfrac-estimate'' uses [[Arithmetic/Rational|rational arithmetic]] and never truncates the intermediate result. When ''terms'' is large, ''cfrac-estimate'' runs slow because numerator and denominator grow big.

math.ranges prettyprint sequences ;

IN: rosetta.cfrac

PRIVATE>

 

{{out}}

<pre> Square root of 2: 1.414213562373095048801688724209

 

=={{header|Felix}}==

let a = match n with | 0 => 3.0 | _ => 6.0 endmatch in

let b = pow(2.0 * n.double - 1.0, 2.0) in

println$ cf_iter 200 sqrt_2; // => 1.41421

println$ cf_iter 200 napier; // => 2.71818

 

=={{header|Forth}}==

{{trans|D}}

: fnapier dup dup 1 > if 1- else drop 1 then s>f dup 1 < if drop 2 then s>f ;

: fpi dup 2* 1- dup * s>f 0> if 6 else 3 then s>f ;

do i over execute frot f+ f/ -1 +loop

0 swap execute fnip f+ \ calcucate for 0

{{out}}

<pre>

 

=={{header|Fortran}}==

implicit none

end function

 

{{out}}

<pre>

 

=={{header|FreeBASIC}}==

 

function sqrt2_a( n as uinteger ) as uinteger

print calc_contfrac( @sqrt2_a, @sqrt2_b, MAX )

print calc_contfrac( @napi_a, @napi_b, MAX )

 

=={{header|Fōrmulæ}}==

 

 

 

 

=={{header|Go}}==

 

import "fmt"

fmt.Println("nap: ", cfNap(20).real())

fmt.Println("pi: ", cfPi(20).real())

{{out}}

<pre>

=={{header|Groovy}}==

{{trans|Java}}

 

import static java.lang.Math.pow

System.out.println(calc(f, 200))

}

{{out}}

<pre>1.4142135623730951

 

=={{header|Haskell}}==

import Data.Char (intToDigit)

 

(putStrLn .

take 200 . decString . (approxCF 950 :: [(Integer, Integer)] -> Rational))

{{out}}

<pre>

3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386

</pre>

 

-- ignoring the task-given pi sequence: sucky convergence

cf2rat_p p s = f $ map ((\i -> (cf2rat i s, cf2rat (1 + i) s)) . (2 ^)) [0 ..]

where

 

-- returns a decimal string of n digits after the dot; all digits should

 

main :: IO ()

 

=={{header|J}}==

 

sqrt2=: cfrac 1 1,200$2 1x

1.4142135623730950488016887242096980785696718753769480731766797379907324784621205551109457595775322165

3.1415924109

 

Note that there are two kinds of continued fractions. The kind here where we alternate between '''a''' and '''b''' values, but in some other tasks '''b''' is always 1 (and not included in the list we use to represent the continued fraction). The other kind is evaluated in J using <code>(+%)/</code> instead of <code>+`%/</code>.

{{trans|D}}

{{works with|Java|8}}

import java.util.*;

import java.util.function.Function;

System.out.println(calc(f, 200));

}

<pre>1.4142135623730951

2.7182818284590455

computes the continued fraction until the difference in approximations is less than or equal to delta,

which may be 0, as previously noted.

def continued_fraction( first; next; count ):

def cf:

if .[0] == count then 0

end;

[2,null] | cf;

'''Examples''':

 

The convergence for pi is slow so we select delta = 1e-12 in that case.

"2|sqrt : \(2|sqrt) : \(continued_fraction_delta( [1,1,1]; [.[0]+1, 2, 1]; 0))",

"1|exp : \(1|exp) : \(2 + (1 / (continued_fraction_delta( [1,1,1]; [.[0]+1, .[1]+1, .[2]+1]; 0))))",

"pi : \(1|atan * 4) : \(continued_fraction_delta( [1,3,1]; .[0]+1 | [., 6, ((2*. - 1) | (.*.))]; 1e-12)) (1e-12)"

{{Out}}

Value : Direct : Continued Fraction

2|sqrt : 1.4142135623730951 : 1.4142135623730951

1|exp : 2.718281828459045 : 2.7182818284590455

 

=={{header|Julia}}==

 

end

end

 

 

{{out}}

 

=={{header|Klong}}==

cf::{[f g i];f::x;g::y;i::z;

f(0)+z{i::i-1;g(i+1)%f(i+1)+x}:*0}

cf({:[0=x;2;x]};{:[x>1;x-1;x]};1000)

cf({:[0=x;3;6]};{((2*x)-1)^2};1000)

{{out}}

<pre>

=={{header|Kotlin}}==

{{trans|D}}

 

typealias Func = (Int) -> IntArray

)

for (pair in pList) println("${pair.first} = ${calc(pair.second, 200)}")

 

{{out}}

pi = 3.141592622804847

</pre>

 

=={{header|Lua}}==

{{trans|C}}

local a = 0.0

local b = 0.0

end

 

{{out}}

<pre>1.4142135623731

 

=={{header|Maple}}==

contfrac:=n->evalf(Value(NumberTheory:-ContinuedFraction(n)));

contfrac(2^(0.5));

contfrac(Pi);

contfrac(exp(1));

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

napier=Function[n,{2,Transpose@{Range[n],Prepend[Range[n-1],1]}}];

pi=Function[n,{3,Transpose@{Array[6&,n],Array[(2#-1)^2&,n]}}];

approx=Function[l,

N[Divide@@First@Fold[{{#2.#[[;;,1]],#2.#[[;;,2]]},#[[1]]}&,{{l[[2,1,1]]l[[1]]+l[[2,1,2]],l[[2,1,1]]},{l[[1]],1}},l[[2,2;;]]],10]];

{{out}}

<pre>

 

=={{header|Maxima}}==

for i from n step -1 thru 2 do z: b[i]/(a[i] + z), a[1] + z)$

 

fpprec: 20$

x: cfeval(cf_pi(10000))$

 

=={{header|NetRexx}}==

* Derived from REXX ... Derived from PL/I with a little "massage"

* SQRT2= 1.41421356237309505 <- PL/I Result

end

Get_Coeffs(form,0)

Who could help me make a,b,sqrt2,napier,pi global (public) variables?

This would simplify the solution:-)

 

=={{header|Nim}}==

var a, b, temp = 0.0

for i in countdown(n, 1):

(a, b)

 

{{out}}

<pre>1.414213562373095

 

=={{header|OCaml}}==

and nap = 2, fun n -> (max 1 (n-1), n)

and root2 = 1, fun n -> (1, 2) in

Printf.printf "sqrt(2)\t= %.15f\n" (eval root2 1000);

Printf.printf "e\t= %.15f\n" (eval nap 1000);

Output (inaccurate due to too few terms):

<pre>sqrt(2) = 1.414213562373095

=={{header|PARI/GP}}==

Partial solution for simple continued fractions.

 

Output:

<pre>%1 = 1.4142135623730950488016887242096980786</pre>

 

=={{header|Perl}}==

 

 

}

 

{{out}}

<pre>√2 ≈ 1.414213562

 

=={{header|Phix}}==

{{Out}}

<pre>

Pi: 3.141592654

</pre>

 

=={{header|PicoLisp}}==

(de fsqrt2 (N A)

(default A 1)

(prinl (format (fsqrt2 200) *Scl))

(prinl (format (napier 200) *Scl))

{{out}}

<pre>

 

=={{header|PL/I}}==

test: proc options (main);

declare n fixed;

put (1 + 1/denom(2));

 

{{out}}

<pre> 1.41421356237309505E+0000 </pre>

Version for NAPIER:

declare n fixed;

 

put (2 + 1/denom(0));

 

<pre> 2.71828182845904524E+0000 </pre>

Version for SQRT2, NAPIER, PI

 

continued_fractions: /* 6 Sept. 2012 */

put skip edit ('PI=', calc(pi, 99999)) (a(10), f(20,17));

 

{{out}}

<pre>

 

=={{header|Prolog}}==

% square root 2

continued_fraction(200, sqrt_2_ab, V1),

pi_ab(0, 3, _).

pi_ab(N, 6, V) :-

{{out}}

<pre> ?- continued_fraction.

=={{header|Python}}==

{{works with|Python|2.6+ and 3.x}}

import itertools

try: zip = itertools.izip

cf = CF(Pi_a(), Pi_b(), 950)

print(pRes(cf, 10))

===Fast iterative version===

{{trans|D}}

 

def calc(fun, n):

print calc(fsqrt2, 200)

print calc(fnapier, 200)

{{out}}

<pre>1.4142135623730950488016887242096980785696718753770

2.7182818284590452353602874713526624977572470937000

3.1415926839198062649342019294083175420335002640134</pre>

 

=={{header|Racket}}==

===Using Doubles===

This version uses standard double precision floating point numbers:

#lang racket

(define (calc cf n)

(calc cf-napier 200)

(calc cf-pi 200)

Output:

1.4142135623730951

2.7182818284590455

3.1415926839198063

 

===Version - Using Doubles===

This versions uses big floats (arbitrary precision floating point):

#lang racket

(require math)

(calc cf-napier 200)

(calc cf-pi 200)

Output:

(bf #e1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358960036439214262599769155193770031712304888324413327207659690547583107739957489062466508437105234564161085482146113860092820802430986649987683947729823677905101453725898480737256099166805538057375451207262441039818826744940289448489312217214883459060818483750848688583833366310472320771259749181255428309841375829513581694269249380272698662595131575038315461736928338289219865139248048189188905788104310928762952913687232022557677738108337499350045588767581063729)

(bf #e2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723624212700454495421842219077173525899689811474120614457405772696521446961165559468253835854362096088934714907384964847142748311021268578658461064714894910680584249490719358138073078291397044213736982988247857479512745588762993966446075)

(bf #e3.14159268391980626493420192940831754203350026401337226640663040854412059241988978103217808449508253393479795573626200366332733859609651462659489470805432281782785922056335606047700127154963266242144951481397480765182268219697420028007903565511884267297358842935537138583640066772149177226656227031792115896439889412205871076985598822285367358003457939603015797225018209619662200081521930463480571130673429337524564941105654923909951299948539893933654293161126559643573974163405197696633200469475250152247413175932572922175467223988860975105100904322239324381097207835036465269418118204894206705789759765527734394105147)

 

=={{header|Raku}}==

(formerly Perl 6)

{{Works with|rakudo|2015-10-31}}

{

my $x = @a[$n - 1];

printf "√2 ≈%.9f\n", continued-fraction(:a(1, |(2 xx *)), :b(Nil, |(1 xx *)));

printf "e ≈%.9f\n", continued-fraction(:a(2, |(1 .. *)), :b(Nil, 1, |(1 .. *)));

{{out}}

<pre>√2 ≈ 1.414213562

 

Raku has a builtin composition operator. We can use it with the triangular reduction metaoperator, and evaluate each resulting function at infinity (any value would do actually, but infinite makes it consistent with this particular task).

map { .(Inf) }, [\o] map { @a[$_] + @b[$_] / * }, ^Inf

}

printf "√2 ≈ %.9f\n", continued-fraction((1, |(2 xx *)), (1 xx *))[10];

printf "e ≈ %.9f\n", continued-fraction((2, |(1 .. *)), (1, |(1 .. *)))[10];

{{out}}

<pre>√2 ≈ 1.414213552

 

More code is used for nicely formatting the output than the continued fraction calculation.

parse arg terms digs .

if terms=='' | terms=="," then terms=500

say right(?,8) "=" $ ' α terms='aT ...

if b.1\==1 then say right("",12+digits()) ' ß terms='bT ...

'''output'''

<pre>

 

===version 2 derived from [[#PL/I|PL/I]]===

* Derived from PL/I with a little "massage"

* SQRT2= 1.41421356237309505 <- PL/I Result

end

call Get_Coeffs form, 0

 

===version 3 better approximation===

* The task description specifies a continued fraction for pi

* that gives a reasonable approximation.

end

call Get_Coeffs 0

 

=={{header|Ring}}==

# Project : Continued fraction

 

eval("temp3=" + expr + "1" + str)

return a0 + b1 / temp3

Output:

<pre>

e = 2.718281828459046

PI = 3.141592653588017

</pre>

 

=={{header|Ruby}}==

 

# square root of 2

puts estimate(sqrt2, 50).to_s('F')

puts estimate(napier, 50).to_s('F')

{{out}}

<pre>$ ruby cfrac.rb

 

=={{header|Rust}}==

use std::iter;

 

println!("{}", continued_fraction!(pia, pib));

}

 

{{out}}

{{works with|Scala|2.9.1}}

Note that Scala-BigDecimal provides a precision of 34 digits. Therefore we take a limitation of 32 digits to avoiding rounding problems.

import Stream._

val sqrt2 = 1 #:: from(2,0) zip from(1,0)

println()

}

{{out}}

<pre>sqrt2:

precision: 3.14159265358</pre>

For higher accuracy of pi we have to take more iterations. Unfortunately the foldRight function in calc isn't tail recursiv - therefore a stack overflow exception will be thrown for higher numbers of iteration, thus we have to implement an iterative way for calculation:

import Stream._

val sqrt2 = 1 #:: from(2,0) zip from(1,0)

println()

}

{{out}}

<pre>sqrt2:

The following code relies on a library implementing SRFI 41 (lazy streams). Most Scheme interpreters include an implementation.

 

(import (rnrs base (6))

(srfi :41 streams))

(stream-cons 3

(stream-cycle (build-stream (lambda (n) (expt (- (* 2 (+ n 1)) 1) 2)))

 

Test:

1.4142135623730951

> (cf->real napier)

2.7182818284590455

> (cf->real pi)

 

=={{header|Sidef}}==

f(func (r) {

r < n ? (a(r) / (b(r) + __FUNC__(r+1))) : 0

for k in (params.keys.sort) {

printf("%2s ≈ %s\n", k, continued_fraction(params{k}...))

{{out}}

<pre>

{{trans|Rust}}

 

@inlinable

public func power(_ n: Self) -> Self {

 

print("π ≈ \(continuedFraction(piA, piB))")

 

{{out}}

e ≈ 2.7182818284590455

π ≈ 3.141592653339042</pre>

 

=={{header|Tcl}}==

{{trans|Python}}

Note that Tcl does not provide arbitrary precision floating point numbers by default, so all result computations are done with IEEE <code>double</code>s.

 

# Term generators; yield list of pairs

puts [cf r2]

puts [cf e]

{{out}}

<pre>1.4142135623730951

 

=={{header|VBA}}==

Private Function continued_fraction(steps As Integer, rid_a As String, rid_b As String) As Double

Dim res As Double

Debug.Print "Napier:", continued_fraction(precision, "nap_a", "nap_b")

Debug.Print "Pi:", continued_fraction(precision, "pi_a", "pi_b")

<pre>Precision: 10000

Sqr(2): 1,4142135623731

=={{header|Visual Basic .NET}}==

{{trans|C#}}

Function Calc(f As Func(Of Integer, Integer()), n As Integer) As Double

Dim temp = 0.0

End Sub

 

{{out}}

<pre>1.4142135623731

2.71828182845905

3.14159262280485</pre>

 

=={{header|XPL0}}==

The number of iterations (N) needed to get the 13 digits of accuracy was determined by experiment.

int N;

real A, B, F;

RlOut(0, 3.0+F); CrLf(0);

RlOut(0, ACos(-1.0)); CrLf(0);

{{out}}

<pre>

 

=={{header|zkl}}==

a0 + [n..1,-1].reduce(

'wrap(p,n){ fb(n)/(fa(n)+p) },0.0) }.fp(fa,fb,a0)

cf creates a function that calculates the continued fraction from the bottom up. The new function takes a single parameter, n, which is used to calculate the nth term.

sqrt2(200) : "%.20e".fmt(_).println();

nap:=cf((0.0).create,fcn(n){ (n==1) and 1.0 or (n-1).toFloat() },2.0);

println(nap(15) - (1.0).e);

pi:=cf((6.0).noop,fcn(n){ n=2*n-1; (n*n).toFloat() },3.0);

(1.0).create(n) --> n, (1.0).noop(n) --> 1.0

{{out}}

=={{header|ZX Spectrum Basic}}==

{{trans|BBC_BASIC}}

20 LET a0=2: LET b1=1: LET a$="N": LET b$="N": PRINT "e = ";: GO SUB 1000

30 LET a0=3: LET b1=1: LET a$="6": LET b$="(2*N+1)^2": PRINT "PI = ";: GO SUB 1000

1035 FOR i=1 TO n: LET p$=p$+")": NEXT i

1040 PRINT a0+b1/VAL (e$+"1"+p$)