I know that given two complex numbers $z_1 = a + bi$ and $z_2 = c + di$, the multiplication of these two numbers is defined as
$$ z_1*z_2 = (ac - bd) + i(ad + cb) $$
I also know that I can easily derive this formula by applying the distributive property of multiplication and considering $i^2 = -1$.
But $i^2 = -1$ is a consequence of the definition of multiplication.
So, the question is: How mathematicians defined the multiplication of complex numbers whether they didn't know the value of $i^2$?
Indeed, $i^2=-1$ is coherent with the definition of multiplication, but things appeared in the reverse order.
While working on the resolution of algebraic equations (roots of polynomials), the need arose to consider square roots of negative numbers. Hence the "imaginary" symbol $i=\sqrt{-1}$, which is enough to represent the square root of any negative number $\sqrt{-a}$ as $\sqrt ai$.
From there, the algebra of complex numbers easily follows, treating $i$ as a variable and using the rule $i^2=-1$:
$$(a+bi)+(c+di)=a+bi+c+di=(a+b)+(c+d)i,$$ $$(a+bi)(c+di)=ac+adi+bci+bdi^2=ac+adi+bci-bd=(ac-bd)+(ad+bc)i.$$
Division is not much more difficult to invent, using a trick to turn the denominator to a real number,
$$\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(a+bi)(c-di)}{c^2+d^2}.$$
Now, what about the square root $\sqrt{a+bi}$ ?
Well, let us try and solve $a+bi=(c+di)^2=c^2-d^2+2cdi$, or $$a=c^2-d^2\\b=2cd.$$ Then, $$a^2+b^2=c^4-2c^2d^2+d^4+4c^2d^2=c^4+2c^2d^2+d^4=(c^2+d^2)^2,$$ so that $$c^2=\frac{\sqrt{a^2+b^2}+a}2,\\d^2=\frac{\sqrt{a^2+b^2}-a}2.$$
As you can see, $c$ and $d$ are always real numbers, so that no new symbol needs to be introduced to take square roots. From the last formulas, we have $$\sqrt i=\frac{1+i}{\sqrt2}.$$
Being fairly new to the site I'm not sure if I can reference a book but there is an excellent text called Visual Complex Analysis by Tristan Needham which explains the origin of complex numbers as a leap of faith required to make a general solution to the cubic equation work. Early mathematicians found that if only they could suppose the existence of this imaginary square root of -1 then they could happily cancel some terms and give the correct solution. It is not as many people suppose to give the equation
$$ x^2 + 1 = 0 $$
a solution because mathematicians could see that it didn't need one. Yves' answer is more to the point but I heartily recommend you check out this book for understanding and fun.
I think the best answer to your question is to distinguish between the historical order in which information developed and the logical order in which that information is now presented.
In the case of complex numbers, the historical order is that people imagined a square root of $-1$, gave it a name $i$, and proceeded to do arithmetic and algebra (and even calculus) with it, using the same computational rules that they already knew for real numbers. In particular, those rules led to the general formula for multiplication that you quoted in the question.
In a rigorous modern presentation of the complex number system, one sets out, right at the beginning, the definitions of the arithmetical operations, including the general formula for multiplication. Then the fact that $i^2=-1$ emerges as just a simple special case of the general definition.
The distinction between the order of discovery and the order of presentation is not a special situation for the complex numbers. For example, in plane geometry, many facts were known before Euclid systematized the subject and deduced vast numbers of facts from just a few axioms. The axioms precede the theorems in a logical development, but that does not mean that the axioms were discovered (and explicitly stated) before the theorems were.
To solve equations, we came up with equations like $x - 1 = 0$, what is $x$?, and the naturals solved them easily. Now, what about $x + 1 = 0$? So we need negative numbers. Wait, what about $2x = 1$? So we need rational numbers. Finally, what about $x^2 = 2$? So we need irrational numbers.
Finally, what about $x^2 = -1$? So we need to invent the "imaginary" numbers to solve it.
Pairs of numbers. Now, we define a new number system, only now there are always pairs of numbers. Then go through explaining how addition and multiplication work.
Then we show that in this system, we have $(0,1)\times (0,1) = (-1,0)$, i.e., we have defined a new system, under which it makes sense to say that $\sqrt{-1} = i$, when $i=(0,1)$. And that's really all there is to imaginary numbers: a definition of a new number system, which makes sense to use in most places. And under that system, there is an answer to $\sqrt{-1}$.
So we have the following definitions: $$(a, b) + (c, d) = (a+c, b+d)$$ $$(a, b)\times (c, d) = (ac-bd, ad+cb)$$
For instance, the programming language Python implements the complex numbers and their algebraic operations. Maybe you can see a good reference in its source code.
By definition $i^2 = -1$. Your concern is about how can one even speak of "$i*i$" at all without having defined the operator "$*$". The answer is: We don't care for now what "$*$" is, as long as there will be a "$*$" product that satisfies $i*i = -1$ and generalizes the notion of product in $\mathbb R$, i.e., $a*c = a \cdot c$, where $a,c \in \mathbb R$ and "$\cdot$" is the real product.
Now you define what is a complex number, and that you can represent them as $z = a + i b$. On top of that some brilliant person defined once the product of two of these numbers to follow also the commutative, associative and distributive laws: $z_1*z_2 = a*c + i(a*d + c*b) + (i*b)*(i*d) = a \cdot c + i(a \cdot d + c \cdot b) + (i*i)*(b \cdot d)$.
Since we defined $i^2 = i*i = -1$, we obtain your result: $z_1*z_2 = (ac - bd) + i(ad + cb)$, as you already know. So, to answer your question, $i^2 = -1$ does not follow from the product rule, but it precedes the notion of the product in $\mathbb C$, we established a definition of "$*$" that is consistent with it.
Everybody agreed on this definition because of the "nice"$^{1}$ properties, e.g. $\mathbb C$ is a ring that generalizes $\mathbb R$, you can represent complex numbers geometrically by vectors in $\mathbb R^2$ and, the nicest property for me, you can divide these vectors by each other!
$^{1}$: Let me elaborate on what is meant by "nice". Let me first give you an example of an "ugly" complex product: Define an operator $**$ such that $z_1 ** z_2 = ac-bd$ and call it "complex quasi-product", it also fulfills $i ** i = -1$ and, moreover, it is consistent with the product of two real numbers!: $z_1 * z_2 = a c$. No one can stop you from using it as a complex multiplication as well, but it lacks of many "nice" (geometrical or algebraic) properties, as, for example, the quasi-product of a real times a complex number won't be to "stretch" the complex vector (in the $\mathbb R ^2$ representation).
I came up with my own derivation of multiplication in the complex field. Maybe it will gain us insight into how earlier mathematicians viewed this subject.
We start with the following question: Find $x$ such that $x^2 +1=0$. Now there does not exist such $x\in\mathbb{R}$ since $x^2 >0,\forall x\in\mathbb{R}$ is an immediate consequence of the Field axioms. Perhaps considering different fields might shine some light on this problem. Let's consider $\mathbb{R}^2$. We try and define addition and multiplication to satisfy the Field Axioms.
Choose $(a,b)\in\mathbb{R}^2$ and $(c,d)\in\mathbb{R}^2$. Then $(a,b)+(c,d)=(a+c,b+d)$ satisfies all the addition axioms. Also, we have that the additive identity element is $\left(0,0\right)$ as a consequence of that definition.
Now, we try and define multiplication. First, we start with defining it as $(a,b)(c,d)=(ac,bd)$. However, $(3,0)\ne(0,0)$ does not have a multiplicative inverse. Therefore, that particular definition breaks the multiplicative inverse Field axiom. It is now clear that any sort of definition of multiplication must involve addition and integrate both elements otherwise elements of $\mathbb{R}^2$ such as $(3,0)$ will always not have an inverse.
We now try the following definition. Define $(a,b)(c,d)=(ac+bd,ad+bc)$. Immediately, we find that the multiplicative identity is $(1,0)$ as a result of that definition. So far good, let's try and find a general equation for the multiplicative inverse of $(a,b)$. We must have that $ac+bd=1$ and $ad+bc=0$. We then end up with $c=\dfrac{a}{a^2-b^2}$ and $d=\dfrac{-b}{a^2-b^2}$. However, this can still break because we are not guaranteed that $a^2 -b^2\ne0$.
We finally end up with the following definition. Define $(a,b)(c,d)=(ac-bd,ad+bc)$. This satisfies all the multiplication Field axioms.
This Field we defined has a solution for $x^2 +1=0$. We have that $(0,1)^2+(1,0)=(-1,0)+(1,0)=(0,0)$. We call this new Field $\mathbb{C}=\mathbb{R}^2$.
