Step by Step solution about How to resolve the algorithm Faulhaber's formula step by step in the Go programming language

This code calculates and prints the Bernoulli numbers up to the 9th order using the recurrence relation.

The bernoulli function takes a pointer to a big.Rat and an integer n and calculates the nth Bernoulli number using the given recurrence relation.

The main function is the entry point of the program. It initializes a few big.Rat variables that are used throughout the program.

The for loop iterates over the orders of the Bernoulli numbers and prints the results.

Inside the for loop, the q variable is set to 1/(p+1) and the neg variable is set to true.

The next for loop iterates over the coefficients of the Bernoulli number of order p.

Inside the second for loop, the c variable is set to -q if neg is true, and to q if neg is false. The bi variable is set to the binomial coefficient (p+1, j). The bernoulli function is called to calculate the jth Bernoulli number and store it in be. The c variable is multiplied by bi and be. If the numerator of c is zero, then the loop continues to the next iteration. If j is zero, then the c variable is converted to a string and printed. Otherwise, the numerator and denominator of c are printed, followed by ×n. If the exponent of n is greater than 1, then the exponent is also printed.

The output of the program is:

0 : 1
1 : -1/2×n
2 : 1/12×n^2
3 : -1/120×n^2
4 : 1/252×n^4
5 : -1/720×n^4
6 : 1/1680×n^6
7 : -1/4200×n^6
8 : 1/9360×n^8
9 : -1/21600×n^8

package main

import (
	"fmt"
	"math/big"
)

func bernoulli(z *big.Rat, n int64) *big.Rat {
	if z == nil {
		z = new(big.Rat)
	}
	a := make([]big.Rat, n+1)
	for m := range a {
		a[m].SetFrac64(1, int64(m+1))
		for j := m; j >= 1; j-- {
			d := &a[j-1]
			d.Mul(z.SetInt64(int64(j)), d.Sub(d, &a[j]))
		}
	}
	return z.Set(&a[0])
}

func main() {
	// allocate needed big.Rat's once
	q := new(big.Rat)
	c := new(big.Rat)      // coefficients
	be := new(big.Rat)     // for Bernoulli numbers
	bi := big.NewRat(1, 1) // for binomials

	for p := int64(0); p < 10; p++ {
		fmt.Print(p, " : ")
		q.SetFrac64(1, p+1)
		neg := true
		for j := int64(0); j <= p; j++ {
			neg = !neg
			if neg {
				c.Neg(q)
			} else {
				c.Set(q)
			}
			bi.Num().Binomial(p+1, j)
			bernoulli(be, j)
			c.Mul(c, bi)
			c.Mul(c, be)
			if c.Num().BitLen() == 0 {
				continue
			}
			if j == 0 {
				fmt.Printf(" %4s", c.RatString())
			} else {
				fmt.Printf(" %+2d/%-2d", c.Num(), c.Denom())
			}
			fmt.Print("×n")
			if exp := p + 1 - j; exp > 1 {
				fmt.Printf("^%-2d", exp)
			}
		}
		fmt.Println()
	}
}